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Research activity carried out by each research associate on the 
Grant NCC2-420 are listed below. The activity covers the period October 
1991 through March 1992. Dr. Greg Wilson is a new associate in this 
grant . 


Dr. Dikran Babikian 

Predicting and characterizing the flowfield of the 20 MW Arc-Jet 
is in progress. Various available computer codes are used for the 
different sections of the flow for initial predict ions . The predictions 
of the free stream spectra has been obtained and compared with the 
experimental data. The prediction of the shock layer spectra and it 1 s 
comparison with the available experimental data is in progress. 


Dx* Tahir Gokcen 

The work in progress is in two areas. First, the coupling of the 
radiative transfer to flow motion for a non-gray emitting and absorbing 
gas has been formulated for two-dimensional nonequilibrium flows. The 
preliminary results will be presented in AIAA 23rd Plasma Dynamics and 
Lasers Conference, July 6-8, 1992, Nashville, Tenn esse. Please see the 



attached copies of the AIAA paper abstract. Second, in collaboration 


with Dr. Iain Boyd of Eloret Institute, thermo chemical models for 
continuum and particle simulations of hypersonic flow is being 
investigated. The results will be presented in AIAA 27th Thermophysics 
Conference, July 6-8, 1992, Nashville, Tennesse. Please see the attached 
copies of the AIAA paper abstract under Appendix A. 


Dr. Seung-Ho Lee 

In previous period, I performed the grid generation for the flow 
field calculations for the NASA-Ames NASP nozzle/afterbody model, which 
had been tested in the Ames 3.5 ft. hypersonic wind tunnel. The result 
of the flow field calculation along with the experimental result 
presented as AIAA Paper 92-0387 in the AIAA 30th Aerospace Meeting & 
Exhibit, Reno, Nevada, January 6-9, 1992, entitled "Single Expansion 
Ramp Nozzle Simulations." 

Pulse facilities , such as the reflected shock tunnel, expansion 
tubes, and gun tunnels, can produce a pulse of high energy airflow over 
a short time, typically on the order of a millisecond. It is this 
ability to produce high energy levels in a test gas that has led to the 
application of these type facilities to the experimental study of 
hypersonic mixing and combustion in a scramjet propulsion system for a 


NASP-type evhicle. Those facilities had been evaluated for the test 
conditions using the various flow codes {mostly one-dimensional codes) 
with their own chemistry models. It is important to evaluate those 
chemistry models and their reactions and choose the best suitable 
model. Application of this model would enable to simulate the 
conditions at the test sections of the pulse hypersonic facilities. 

The result of this study will be presented in the AIAA 17th Aerospace 
Ground Testing Conference, Nashville, Tennessee, July 6-8, 1992, 
entitled "Thermo chemical Relaxation in High Enthalpy Nozzle Flows” (AIAA 
Paper 92-4015) . Please see the attached copies of the AIAA paper 
abstract under Appendix B. 


Dr. Dinesh Prabhu 

The following tasks were initiated and completed for the reporting 
period Oct . 1 91-Mar . 1 92 - (1) As a preliminary test of the capability of 
the computer codes, the flow past a two-dimensional wedge was simulated 
using nonequilibrium upwind codes developed by Dr.J.-L. Cambier. The 
computed results were compared with those obtained at Langley Research 
Center. The good agreement between the codes gave the confidence to 
compute the flow in the inlet of the NASP combustor. (2) The performance 


of the NASP combustor inlet (modelled as a two-dimensional geometry) was 


studied for the 16" Shock Tunnel flow conditions corresponding to driver 
pressures of 6000 psi and 8000 psi. The effect on the performance due to 
varying cowl lengths and flow Mach numbers was also investigated. The 
location of the impingement point of the cowl shock was determined for 
various cowl lengths and special attention was paid to the case of shock 
cancellation at the body surface. This work was done in collaboration 
with Dr. J.-L. Cambier whose two-dimensional/axisymmetric nonequilibrium 
codes were used and during the course of this study a detailed grid 
sensitivity analysis of the computer codes was done. The adaptive gird 
code SAGE developed by Dr .Venkat apathy and Dr. Davies was also used in 
order to more accurately obtain the impingement points. (3) In support 
of the absorption experiment to be located at station N3 of the 16" 
shock tunnel nozzle, the NO (y) system [(0,0) and (0,1) bands] was 
studied in absorption with both quasi-one dimensional and two- 
dimensional axisymmetric reacting flow analyses. The required depth of 
the optical access, to minimize the effects of the hot boundary layer on 
the measurements, was estimated. The NO (y) band system was studied in 
very high-resolution with and without splitting of the 2 Z- 2 n 
transitions. The entire analysis was done with the assumption of 
chemical nonequilibrium and thermal equilibrium using a spectral code 
NEQAIR developed by Dr. Ellis Whiting. (4) A similar analysis was done 
at the exit of the nozzle in order to characterize the chemical state of 
the gas. (5) The spectra code was modified in order to do the 
spectroscopic study of the contaminants in the shock tunnel. The 



spectroscopic data for the species Fe, Cr, Ni and Cu are now 
incorporated into the code and the code has been handed over to the 
experimentalists who plan on using it to study the contamination levels 
in the shock tunnel. (6) The calibration of the Pinckney static pressure 
probes was initiated using a nonequilibrium parabolized Navier-Stokes 
code. Please see the Appendix C for further details. 


Ms* Susan Tokarcik 

1. I participated in the Ames Technical Paper Contest for Women. 

I presented work on hypersonic drag devices . I was one of the winners 
of this contest and will be attending the Society of Women Engineers 
Conference in June to present my winning paper. Also, a paper on this 
material has been submitted to the Journal of Spacecraft and Rockets for 
publication . 

2. I have also been working on computing the flowfield around a 
lunar return aerobrake concept proposed by Mike Tauber. I am working 
with E. Venketapathy on this project. Thus far an axisymmetric study of 
the aerobrake design has been completed for an aerobrake with flares at 
angles of 50°°, 55°°, and 60°°. Three-dimensional calculations have also 
been completed at the same flare angles. It has been found that because 
of the large separated region created by the flares, the flow about the 
aerobrake is very complicated. Further investigation of this flow is 


anticipated and will include calculation of the base flow region. 
Comparisons between the axisymmetric solutions and the 3-D solutions 
have been made for surface pressure, size of the separated region, and 
total drag produced. Validation of the drag calculations and general 
flowfield character are to be completed by calculating a similar 
flowfield that has been studied experimentally. A representative 
picture of this work is shown in Figure 1. 

3. I am continuing as point of contact for the University of 
Maryland. Maryland is involved in the CFD support effort for the Ames 
16” shock tunnel. Thus far J.L. Cambier's axisymmetric TVD code 
(Emozart_2dA. f ) has been transfered to the University of Maryland as 
well as the configuration of the ARC 16” shock tunnel nozzle and a 
generic combustor design. Cambier's 3-D code will eventually be 
transfered to the University of Maryland so that they can perform a full 
3-D calculation of the test section flow with the combustor model in the 
test section. This study is being performed in order to determine the 
likelihood of choking the tunnel with the combustor model. 

4. My involvement with the CFD support of the 16" shock tunnel 
has increased in the last six months and will increase further when D. 
Prabhu leaves the project. I am also working with J.L. Cambier on this 
project A fully turbulent, nonequilibrium chemistry calculation of the 
nozzle flow has been completed. Information from this solution was used 



to characterize the flow in the test section and also to determine the 


optical character of the flow in the nozzle. A study of grid effects on 
the start-up solution has also been performed. This study was done 
assuming an inviscid, non-reacting flow. This start-up process begins 
at the breaking of the the diaphragm between the driven tube and the 
nozzle throat. A representative figure demonstrating grid effects is 
shown in Figure 2. A study of the entire start up procedure for a "best 
case" grid was also completed. This was also for inviscid, non-reacting 
flow. Figure 3. shows a time history of the start-up procedure. 

A study of the applicability of calculating quasi-one-dimensional 
(Q-1D) flow of a reacting gas in order to simulate the full 
axisymmetric, turbulent solution was performed. It is shown that by 
decreasing the nozzle radius by 13% (an 87% nozzle) , the flow in the 
turbulent, reacting case could be simulated successfully at a fraction 
of the CPU time required for the axisymmetric calculation. The results 
of this study are shown in Figure 4. Note that the Q-1D solution gives 
no information about boundary layer thickness or core flow profiles. 

This study will be continued further in order to determine if some 
"rules of thumb" exist for simulating axisymmetric, turbulent flows with 
finite rate chemistry with a Q-1D code with finite rate chemistry. The 
Q-1D code is also being used to study the effects of varying the 
chemical reaction rates for the reactions taking place in the nozzle. 



Another study being done at this time is the simulation of the 
full driver tube, driven tube and nozzle flows in a time accurate 
manner. This simulation is axisymmetric, turbulent, and has finite rate 
chemistry. A representative figure of this flow is shown in Figure 5. 

In this figure, the diaphragm between the driver tube and the driven 
tube has partially ruptured and a shock has formed in the driven tube. 

The results completed thus far for the Ames 16" shock tunnel will 
be presented at the AIAA Thermophysics Conference in July. 

Please see the Appendix D for further details. 


Dr. Ethira-i Venkatapathy 

The axis singularity problem for three dimensional blunt body 
flows were addressed and two possible methods were found to work equally 
well. These methods and the results were reported at the 30th AIAA 
Aerospace Sciences Meeting in January 1992. The work was performed in 
collaboration with Mr. Grant Palmer of NASA Ames research Center. A 
copy of the publication entitled, "Effective Treatment of the Singular 
Line Boundary Problem for Three Dimensional Grids, " is attached at the 


end. 



Related to the National Aero-Space Plane (NASP) Program, 
simulations of the single expansion ramp nozzle experiment were 
completed and the results along with preliminary experimental 
comparisons were reported in a AIAA publication entitled " Single 
Expansion Ramp Nozzle Simulations, " at the 30th Aerospace Sciences 
Meeting in Reno in January 1992. This work was done inn collaboration 
with Mr. S.M. Ruffin of NASA Ames, Dr. E. Keener of Eloret Institute, 
Dr. F. Spade of McDonnel Douglas Research Laborataries of St. Louis and 
Dr. Seung-Ho Lee of Eloret Institute. A copy of the publication is 
attached at the end. 


To understand the flow around a Lunar Return Aero-Brake Vehicle 
and to predict total drag for a number of configurations, flow 
simulations were performed with an axisymmetric and full 3-D geometries . 
Preliminary predictions indicate an effective way of predicting drag for 
a number of 3-D configurations with just axisymmetric computations. 
Though flow details for the full 3-D configurations are very complex and 
completely three dimensional, integrated quantities such as drag 
coefficient can be evaluated with simpler computations. These findings 
are very useful in the design process. This work was performed in 


collaboration with Ms. Susan Tokarcik of Eloret Institute and Dr. Mike 



Taubar of NASA Ames Research Center. Complete computations and results 
on this topic will be presented in an upcoming AIAA conference. 

Please see the Appendix E for further details. 


Dr. Greg Wilson 

Work for this reporting period has involved a transition from my 
Ph.D. research at Stanford University to research involving the NASA 
Ames 16 inch, combustion driven shock tunnel. My dissertation was 
published in December and an additional work co-authored by Myles A. 
Sussman, a Stanford University graduate student, and myself was 
submitted to the AIAA Journal. Both works present numerical simulations 
used to investigate periodic combustion instabilities observed in 
ballistic-range experiments of blunt bodies flying at supersonic speeds 
through hydrogen-air mixtures. The computations are validated by 
comparing experimental shadowgraphs to shadowgraphs created from the 
computed flowfields and by comparing the frequency of the experimental 
and computed instabilities. The simulations add insight into the 
physical processes observed in the experiments and help establish the 
usefulness of the ballistic-range experiments for the validation of 


numerical methods and chemical kinetic models. 



Initial research involving the Ames 16 inch Shock tunnel has 
involved time-dependent quasi-one-dimensional simulations using finite 
rate chemistry to simulate the driver and driven sections of the 
facility. Studies to date have helped quantify viscous and heat 
conduction losses during shock tunnel runs and to demonstrate the effect 
of driver tube pressure gradients on tunnel performance. Numerical 
simulations which include these phenomena compare well with experimental 
data obtained from the facility. 

To compliment the numerical simulations of the shock tunnel, 
expansion tube flow is also being examined using the same numerical 
tools. This is in support of work being done by Dr. Seung-Ho Lee. 

Wilson, G. J., ''Computation of Steady and Unsteady Shock-Induced 
Combustion Over Hypervelocity Blunt Bodies, 1 1 Ph.D Thesis, Stanford 
University, December, 1991. 

Wilson, G. J. and Sussman, M. A., ''Computation of Unsteady Shock- 
Induced Combustion Using a Logarithmic Form of the Species Conservation 
Equations, 11 Submitted to the AIAA Journal, December, 1991. 


Please see the Appendix F for further details. 
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Abstract submitted for 23rd AIAA Plasma Dynamics and Lasers Conference 
July 6-8, 1992, Nashville, Tennessee. 


THE COUPLING OF RADIATION TO BLUNT BODY FLOWS 
IN THERMOCHEMICAL NONEQUILIBRIUM r 

by 

Tahir Gokgen f 

Eloret Institute, 3788 Fabian Way, Palo Alto, CA 94303 
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In hypersonic flight, the radiative energy transfer to the surface of a blunt reentry 
body becomes an important mode of heat transfer as very high temperatures are attained 
behind the bow shock. The radiative transfer occurs along with the other modes of heat 
transfer, conduction find convection. For most of the aerothermodynamic flow problems, 
the radiation and flowfield are not coupled, i.e., the effects of radiation emission and ab- 
sorption on flowfield are neglected. The radiative energy flux to the wall is then calculated 
from this uncoupled flowfield solution. On the other hand, for flow problems such as pre- 
dicting aerobrake flow environment of planetary entry vehicles, radiative transfer becomes 
significant enough that its presence alters the flowfield, which in turn affects the radiative 
heat transfer at the wall. In such situations the radiation and flowfield are said to be 
strongly coupled, therefore, a simultaneous treatment is necessary. This coupling leads to 
a set of integro- differential equations, which are much more complicated than the governing 
differential equations of non-radiating flows. 

Main objective of the proposed paper is to present the numerical tools to simulate two 
dimensional or axisymmetric thermo chemical nonequilibrium flows with strong radiation 
coupling. In particular, a thermo chemical nonequilibrium radiating shock layer about a 
blunt body is considered, and a schematic of such flowfield is shown in Figure 1. 

The present nonequilibrium gas model for air consists of eleven chemical species, (iV 2 , 
0 2 , NO , iV, 0, N 2 , e“), and the thermal state of the gas is described 

by three temperatures: translational, rotational and vibrational (vibrational-electronic). 
The radiative energy transport is assumed to be one- dimensional, the dimension being 
normal to the wall. This is the so-called tangent-slab approximation which is applicable 
to the stagnation region of a blunt body or to blunt bodies of large radius of curvature. 
The thermal radiation is assumed to be governed by the vibrational temperature. 1 The 
gases are assumed to be non-gray, and the spectral absorption coefficients of the species 
are prescribed by 615 absorption bands. 2 

f Research Scientist, Member AIAA 
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Fig* 1* Schematic of a nonequilibrium radiating shock layer about a blunt body. 

The governing Navier-Stokes equations are augmented with the equations accounting 
for thermochemical nonequilibrium processes. The equation set consists of sixteen par- 
tial differential equations: eleven mass conservation equations for species, two momentum 
equations for two dimensioned flows, and three energy equations. The numerical approach 
to solve the governing equations follows the previous work on non-reacting and reacting 
flows, and incorporates severed features developed eeirlier. 3-5 The method is fully implicit 
for fluid dynamics, chemistry, and radiative transfer. It uses flux vector splitting for con- 
vective fluxes, and shock capturing edong with adaptive grid strategy is also implemented. 
The radiative flux is treated fully implicitly, which results in the full block matrix differ- 
ence equations. The full block matrix is inverted iteratively using block iteration methods 
such as Gauss-Seidel Line Relaxation. 6 

The development of the code for two dimensional blunt body flows with fully coupled 
radiation is currently in progress. Therefore, preliminary results are presented. As a test 
case, hypersonic ionizing flow over a 45° blunt cone is considered. Freestream conditions 
of Uqo = 12.km/s, Too = 300. °K , and /><» = 1. x 10~ 4 kg/m 3 are prescribed. The selected 
velocity is high enough that both the ionization and the radiation phenomena play impor- 
tant role in the flowfield behind the shock wave. Figure 2 depicts the blunt cone geometry 
and computational adaptive grid. For these computations flow is assumed to be inviscid, 
and the grid is only adapted to the shock wave. As an example of the flowfield chemistry, 
the electron mole fraction contours near the nose region are plotted in Figure 3, showing 
the extent of ionization. The electron mole fraction and temperature profiles along the 
stagnation streamline are presented in Figure 4 and 5, respectively. In Figure 4, there are 
two distinct peaks in the electron mole fraction plot. The first peak is attributed to the 
associative ionization processes and the second to the electron-impact ionization. 7 
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Fig. 2. Computational grid (50 X 30) over a 45° blunt cone. 



X 

Fig. 3. Electron mole fraction contours near the nose region of the 45° cone. Freestream 
conditions: u* = 12.km/s, = 300.° K, = 1. x 10~ 4 kg/m S . 
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Fig. 4. Electron mole fraction profile along the stagnation streamline of the 45° cone. 
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Fig. 5. Temperature profiles along the stagnation streamline of the 45° cone. 

This shift in the dominant ionization process along the streamline in turn affects the 
temperature profiles. As seen from Figure 5, the temperature profiles show an inflection 
point in the relaxation zone rather than monotonically approaching to the equilibrium 
value. Similar behavior is observed by Park’s one dimensioned computations. 7 For radiative 
transfer computations of the test case, the spectral absorption coefficient k\ is assumed 
to be constant (tt\ = 1.0 m -1 , gray gas). The variation of toted radiative flux along the 
stagnation streamline is plotted in Figure 6. The radiative flux is constemt upstream of 
the shock since the upstreeun gas is specified to be transparent. 

As mentioned earlier, these results are preliminary, the proposed paper will give results 
for both gray and non-gray gases in much greater detail. The emphasis will be on the 
numerical solutions for nonequilibrium blunt body flows where the flowfield and radiation 
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field are strongly coupled. The paper will also present details of the numerical procedures 
associated with the radiation coupling. 



X/R 

Fig. 0. Variation of the toted radiative flux along the stagnation streamline of the 45° cone. 


References 

1 Park, C., Nonequilibrium Hypertonic Aerothermodynamics , John Wiley and Sons, Inc., 
New York, 1989. 

2 Park, C., and Milos, F. S., “Computational Equations for Radiating and Ablating 
Shock Layers,” AIAA paper 90-0356, 1990. 

3 MacCormack, R. W., “Current Status of the Numerical Solutions of the Navier-Stokes 
Equations,” AIAA paper 85-0032, 1985. 

4 Candler, G.V.,~“The Computation of Weakly Ionized Hypersonic Flows in Thermo- 
Chemical Nonequilibrium,” Ph. D. Thetis, Stanford University, 1988. 

5 Gokgen, T., “Computation of Hypersonic Low Density Flows with Thermo chemical 
Nonequilibrium,” Ph. D Thesis, Stanford University, 1989. 

6 Gokgen, T., and Park, C., “The Coupling of Radiative Transfer to Quasi 1-D Flows 
with Thermochemical Nonequilibrium,” AIAA paper 91-0570, 1991. 

7 Park, C., Howe, T. H., Jaffe, L. R., and Candler, G. V., “Chemical-Kinetic Problems 
of Future NASA Missions,” AIAA paper 91-0464, 1991. 


5 




Abstract submitted for consideration to 


AIAA 27th Thermophysics Conference, July 6-8 1992 
Nashville, Tennessee 

i 

ASSESSMENT OF THERMOCHEMICAL MODELS FOR 

CONTINUUM AND PARTICLE SIMULATIONS OF HYPERSONIC FLOW 

*• 

. . s - 
Iain D. Boyd and Tahir ,Gok$en 

s'- 

Eloret Institute 
3788 Fabian Way 
Palo Alto, California. 

Computations are presented for strong shock waves that axe typical of those that form 
in front of reentering spacecraft. The fluid mechanics and thermochemistry are modeled us- 
ing two different approaches. The first employs traditional continuum techniques in solving 
the Navier-Stokes equations. The important chemical effect of vibration-dissociation cou- 
pling is simulated using a multi-temperature model. The second approach employs a parti- 
cle simulation technique (the direct simulation Monte Carlo method, DSMC). This method 
employs a vibrationally-favored dissociation probability model. These techniques are as- 
sessed through comparison of solutions generated under hypersonic flow conditions. Sepa- 
rate cases axe considered that are dominated by thermal relaxation processes (translational- 
rotational-vibrational energy exchange), and chemical relaxation processes (dissociation- 
ionization). These phenomena are each considered for pure nitrogen and comparison made 
with limited experimental data. The studies performed provide an improved understand- 
ing of the relationship between the thermochemical models employed in continuum and 
particle simulations of hypersonic flow. 

A space-vehicle passing through the earth’s atmosphere will traverse a number of 
different flow regimes. At lower altitudes, the fluid density is sufficiently large for the flow 
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to be considered as lying in thermochemical equilibrium. However, as the vehicle ascends 
higher into the atmosphere, the molecular collision rate will fall, and low-density effects 
will become increasingly important. 

Continuum methods are applied to flows in which the collision rate -of the gas is 
sufficient to maintain Boltzmann energy distributions for the various thermal modes of the 
gas. It is not necessary that the temperatures associated with each of the different modes 
be equal, or that chemical equilibrium prevails. As -a general, and somewhat arbitrary 
guideline, it is usually assumed that this condition is' met when the Knudsen number of 

t ■ • \ .♦ 

the flow lies below about 0.01. For computation in the transition flow regime, which lies 
beyond this continuum limit, the most generally successful approach has been the particle 
simulation technique such as DSMC. 

The computation of flow properties for the flight trajectories of many space vehicles 
will require the use of both solution techniques mentioned above. The interface between 
the different flow regimes is therefore of great importance. Clearly, it is desirable to ob- 
tain consistent results with these numerical methods in an overlapping near-continuum 
flow regime. Although the thermochemical models employed in continuum and particle 
methods are quite different., under conditions of thermochemical equilibrium they should 
provide identical solutions. The relationship between a continuum solution and a particle 
simulation under conditions of thermochemical nonequilibrium, however, has not been in- 
vestigated previously. It is therefore the purpose of the present paper to study these effects 
by computing typical hypersonic flows with both continuum and particle simulation meth- 
ods. The continuum and particle approaches employed in this work are briefly described 
below. 

The continuum nonequilibrium gas model for air consists of eleven chemical species, 
(IV 2 , O 2 , NO, N , O, N 2 , O 2 , NO + , N + , 0 + , e - ), and the thermal state of the gas 
is described by three temperatures: translational, rotational and vibrational (vibrational- 
electronic). The governing Navier-Stokes equations are supplemented by the equations 
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accounting for thermochemical nonequilibrium processes. The equation set consists of 
sixteen partial differential equations: eleven mass conservation equations for species, two 
momentum equations for two-dimensional flows, and three energy equations. The numer- 
ical approach to solve the governing equations is described in Refs. 1-3 in. detail. The 
method is fully implicit for fluid dynamics and chemistry. It uses flux vector splitting for 
convective fluxes, and shock capturing. An adaptive grid strategy is also implemented. 
For the computations in this paper, a two-dimensional axisymmetric blunt-body code is. 

>. S 5 t \ 

used and a freestream of pure nitrogen is prescribed... • 
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The particle simulation code employed in this investigation provides modeling of the 
translational, rotational, vibrational, and electron energy distributions. These are comple- 
mented through simulation of dissociative, recombinative, ionizing, and exchange reactions. 
The code is vectorized for efficient execution on a Cray-YMP. Description of the vectorized 
implementation may be found in Refs. 4 and 5. The inclusion of electrons in the simula- 
tion is discussed in detail in Ref. 6. The code employs energy- dependent probabilities of 
rotational and vibrational energy exchange, and a vibrationally-favored dissociation model 
which are all described in Ref. 5. The thermal relaxation rates employed in the particle 
simulation have been equated with the continuum values through the definition provided in 
Ref. 7. The chemical rate data was chosen to correspond to the continuum values wherever 
possible. — 

Computations Rave been performed for three different sets of flow conditions for pure 
nitrogen and these are listed in Table 1. The flow conditions for case 1 correspond to those 
investigated experimentally by Sharma 8 . The density and temperature for cases 2 and 3 
have been chosen to correspond to the ambient conditions of the earth’s atmosphere at an 
altitude of 78 km. The freestream conditions provide increasing enthalpy corresponding 
to flow which is dominated in case 2 by thermal relaxation processes, and in case 3 by 
chemical relaxation processes. 

One of the flow quantities which is most sensitive to changes in the modeling of 
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thermochemistry is the vibrational temperature. This quantity is also of great physical 
significance due to its importance in the determination of radiative emission. Several 
computed profiles of vibrational temperature for case 1 are compared in Fig. 1. The 
continuum profile was obtained with Park’s two- temperature dissociation model in which 
the controling temperature is given by the square root of the product of the translational 
and vibrational values. Three different solutions have been obtained with the particle 
method. The vibrationally-favored dissociation model (VFD) employed in the particle 
simulations computes a dissociation probability based on the total collision energy ’(the 
sum of the translational, rotational, and vibrational components) and may be biased in 
addition to the vibrational energy (see Ref. 5). The particle solution labeled (#1) employed 
the total collision energy to evaluate the dissociation probability, and set the degree of 
additional vibration-dissociation coupling to zero. This is the configuration which has 
been in general usage and leads to a dissociation rate which is higher than that computed 
with the continuum model. Hence the particle computed vibrational temperatures lie 
well below the continuum values. The particle solution labeled (#2) employed the total 
collision energy together with the degree of vibration-dissociation coupling determined 
in Ref. 5 for nitrogen. In this case, the particle simulated dissociation rate is slightly 
slower than the continuum rate. Finally, the particle solution labeled (#3) employed a 
total energy given_by the sum of the translational collision energy, and the rotational 
and vibrational energies of only the dissociating molecule, and the vibration-dissociation 
coupling is retained. This model has the effect of speeding up the dissociation rate, and 
gives a profile for vibrational temperature which is in good agreement with the continuum 
computations. This configuration also represents the most physically realistic methodology 
for simulating dissociation within the framework of a phenomenological model. 

The temperature profiles of the three energy modes computed with the continuum 
technique and particle method, for configuration (#3), are shown in Fig. 2. The ex- 
perimental data of Sharma 8 is also included. Good agreement between the particle and 
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continuum solutions for all three temperature profiles is obtained. The results obtained for 
the rotational and vibrational energy modes offer good correspondence to the experimental 
values both in the shock-front and in the equilibrium region far behind the shock. The 
close accord between the numerical methods is further demonstrated in Fig. 3 where the 
profiles for the mass fraction of atomic nitrogen axe observed to compare well. 

The density profiles computed with the two numerical methods for case 2 are compared 
in Fig. 4. The continuum solution is shown as a'solid lipe, and the particle solution-is 
shown as a dotted line. The thermochemipal model' employed in the particle simulation 
was configuration (#3). The distance through the shock-wave has been normalized by the 
nose radius (RN=2.3 m),. The density rise predicted in the particle simulation preceeds that 
computed by the continuum method, producing a thicker shock-wave. In the relaxation 
zone behind the shock, there is generally good agreement between the two solutions. The 
computed temperature profiles are compared in Fig. 5. The rise in each of the temperatures 
computed with the particle method also preceeds that in the continuum solution. These 
types of differences between particle and continuum solutions have been found previously 
by Lumpkin and Chapman 6 who noted that the continuum equations predict shock-wave 
thicknesses which are smaller than those measured experimentally and which are predicted 
successfully with DSMC. The differences occur because the continuum formulation breaks 
down in regions of" the flow where high gradients give rise to relatively large values for 
the locaTKnudseri number. In the present study, this occurs in the shock-front, where the 
temperatures first begin to rise. 

Comparison is made of the species mole fractions computed with the two numerical 
methods in Fig. 6 for the 10 km/s conditions of case 3. Only the 3 major species are 
represented. Once again, the particle simulation was performed with the thermochemical 
model labeled (#3). Generally, quite good agreement is found for the two solutions. 
However, as in case 2, it is found that the thermochemical relaxation processes computed 
with the particle method preceed the continuum simulation results due to the presence of 
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a thicker shock. 

The good agreement found between the continuum and particle computations for 
case 1 is partly due to the relatively high-density freestream conditions which place the 
flow in the near-continuum regime. For cases 2 and 3, in understanding the differences 
between the continuum and particle solutions, it is difficult to separate the effects of 
chemical nonequilibrium from those of thermal nonequilibrium. Further studies are being 
undertaken to investigate these coupled phenomena in greater detail, and will be reported 
in the final paper. 


Table 1. Flow conditions. 


Case 

Uoo (km/s) 

Poo (kg/m 3 ) 

Too (K) 

1 

6.2 

1.50xl0~ 3 

300 

2 

5 

2.75xl0“ 5 

188 

3 

10 

2.75xl0~ 5 

188 
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THERMOCHEMICAL RELAXATION IN HIGH ENTHALPY NOZZLE FLOWS 

Seung-Ho Lee, Eloret Institute, Palo Alto, CA 
Mark Loomis, MCAT Institute, Moffett Field, CA 
Chul Park, NASA Ames Research Center, Moffett Field, CA 

Summary 

The problem of thermochemical relaxation in a high enthalpy wind tunnel nozzles is re- 
viewed. The operating conditions of the shock tunnels and axc-jet wind tunnels existing 
in the United States and elsewhere are first surveyed. Next the existing one-dimensional 
computer codes that calculate the thermochemical nonequilibrium state in the nozzle are 
reviewed and compared. The codes are run for the experimental conditions of CALSPAN 
and elsewhere in which vibrational and electron temperature and electron density have 
been measured. The discrepancy between the calculation and the experimental data is 
discussed. Comparison is made also between the interferometric fringe pattern obtained in 
a flow over a wedge placed in the 16 inch Combustion-Driven Shock Tunnel at Ames Re- 
search Center and that predicted by a code. A set of reaction rates is selected to reproduce 
the experimental data. Using this validated code, the degree of simulation is calcualted 
for the existing shock tunnels and arc-jet wind tunnels. 

Introduction 

With the advent of the aerospace vehicles such as Space Shuttle, an effort has been ex- 
pended in recent years to produce a high enthalpy flow in a ground-based facility for testing 
of the models for such vehicles. High enthalpy flows are produced either by means of an 
axc-jet wind tunnel, hot-shot wind tunnel, expansion tube, or shock tunnel. To determine 
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the thermochemical state of the flow produced in such facilities, several computer codes 
have been written. These codes require several thermochemical rate parameters as inputs. 
Validity of these codes and the rate parameters used have not yet been made, mostly 
because required measurements of the flow characteristics are difficult. Each laboratory 
uses different computer codes to determine the facility in the laboratory. Thus, there are 
several questions regarding the performance of such facilities. They are: 

(1) What actually happens in the nozzle? 

(2) Which computer code most closely reproduces the flow conditions in the nozzle? 

(3) How do various facilities compare in their performance? 

The purpose of the present work is to answer these questions. 

What Actually Happens in the Nozzle? 

The actual thermochemical phenomena in a high enthalpy nozzle flow can be understood 
by analyzing the experimental data obtained in various shock tunnels. The first of these 
sets of experimental data are those obtained at Cornell Aeronautical Laboratory in the late 
1950s and early 1960s 1 in which vibrational temperature has been measured. The second 
set of such data are those taken also in Cornell Laboratory in which electron temperature 
and electron density have been-measured. 2-4 The third set of data is obtained at Ames 
Research Center in which the electron density, electron temperature, and electronic exci- 
tation temperature have been measured. 5,6 Lastly, very recently, interferometric density 
measurement has been made at Ames Research Center in a shock tunnel for a flowfield 
produced by a two-dimensional wedge for both air and nitrogen. The interferogram, re- 
produced in Fig. 1, contains a peculiar bending of fringes, which indicates that density is 
decreasing toward the wall. This phenomenon is interpreted to be due to atomic recombi- 
nation. The extent of the bending of the fringes will depend on the amount of dissociated 
atoms in the test section. By analyzing this data, it should be possible to deduce the 
degree of dissociation in the test section flow. 
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In the present work, a comprehensive one-dimensional computer code, named NOZ3T, has 
been developed recently to numerically reproduce the above-mentioned experimental data. 
The code allows three temperatures, namely, the translational-rotational, vibrational, and 
electron-electronic, to be different. This was necessitated by the fact that some of the 
experimental data axe in the form of vibrational temperature or electron temperature. An 
example of the use of this code is shown in Fig. 2, in which the experimental data of Dunn 
and Lordi 2 is reproduced. Such comparison is being made presently for all other existing 
experimental data. 

A two-dimensional, thermochemical nonequilibrium flow calculation is currently being un- 
dertaken to reproduce the interferogram of Fig. 1 by using the CENS2H code developed 
by Park and Yoon. 7 The degree of dissociation of the test section flow is being varied to 
find the value which best reproduces the interferogram. The results are not yet available, 
but will be presented in the final paper. The reaction rate set that produces the observed 
degree of dissociation will be determined by running the NOZ3T code. This completes 
validation of the code N0Z3T and the reaction rate set. 

Which Code and What Rate Set Best Reproduce The Flow? 

The validity of the existing computer codes that calculate the relaxation phenomena in a 
high enthalpy nozzle flow and the associate rate parameters is tested by comparing against 
the N0Z3T code so validated. The code so compared axe: (1) CALSPAN code, 8 (2) Lewis 
GCKP code, 9 and (3) CHEMKIN code. 10 The comparison has not yet begun because 
the validation of N0Z3T is not yet been finished. The result of the comparison will be 
presented in the final paper. 

How Do Each Facilities Perform? 

Using the above-mentioned computer codes, comparison is being made between various 
different facilities in the world. The facilities chosen for this comparison are, for the 
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present: (1) Ames 16 inch combustion-driven shock tunnel, (2) the piston-compression- 
heated shock tunnel T4 in Australia, (3) the piston-compression-heated shock tunnel T5 
in California Institute of Technology, (4) the expansion tube at GASL. More facilities will 
be included in this list before the final paper. In Table 1, first of such comparisons is shown 
between the Ames 16 inch shock tunnel and T4. 

Conclusions 

It is the intention of this study to draw conclusion on (1) what computer code should be 
recommended for determining the thermochemical state in a nozzle in the high enthalpy 
regime, and (2) how each facility compare in their performance. 
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Table 1. Comparison of the flow conditions for Ames 16 inch shock tunnel and T4. 



Ames 16” Shock Tunnel 

T4 (Mach 4) 

T4 (Mach 10) 

Stagnation Chamber Condition 




Temperature (K) 

5750 

8000 

8000 

Pressure (mpa) 

24.13 

30.4 

30.4 

Enthalpy (mj/kg) 

9.630 

16.228 

16.228 

Nozzle Exit Condition 




Location from throat (m) 

5.728 

0.512 

1.777 

A/A* 

270.8 

28.0 

624.5 

Temperature (K) 

678.8 

2992.0 

632.5 

Pressure (kpa) 

1.814 

11.348 

0.477 

Density (kg/m 3 ) 

0.00925 

0.01235 

0.00227 

Velocity (m/sec) 

4136.3 

4586.5 

4988.5 

Mach number 

8.13 

4.24 

5.43 

Species Concentrations 




E" 

1.73xl0 -13 

4.96xl0 -12 

8.44xl0" 12 

N 2 

7.27xl0 -1 

7.31xl0 -1 

7.37 xlO -1 

O 2 

1.93xl0 -1 

1.27xl0 -1 

3.56xl0 -2 

Ar 

1.29xl0 -2 

1.29xl0 -2 

1.29xl0 -2 

N 

8.75xl0 -12 

2.52xl0 -5 

4.94xl0" 8 

0 

5.98xl0“ 3 

7.73xl0 -2 

1.75xl0 -1 

NO 

6.19xl0~ 2 

5.13xl0" 2 

3.95xl0 -2 

N0+ 

9.80xl0" 9 

2.72xl0 -7 

4.62xl0 -7 








Figure 1. Holographic interferogram over a wedge at 11° inclination in the Ames 16 inch 
shock tunnel flow at H = 10 MJ/kg; test gas = air, reflected-shock pressure = 400 psia, 
with hydrogen injection. 
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The following tasks were initiated and completed for the reporting 
period Oct . ' 91-Mar . 1 92 - (1) As a preliminary test of the capability of 
the computer codes, the flow past a two-dimensional wedge was simulated 
using nonequilibrium upwind codes developed by Dr.J.-L. Cambier. The 
computed results were compared with those obtained at Langley Research 
Center. The good agreement between the codes gave the confidence to 
compute the flow in the inlet of the NASP combustor. (2) The performance 
of the NASP combustor inlet (modelled as a two-dimensional geometry) was 
studied for the 16" Shock Tunnel flow conditions corresponding to driver 
pressures of 6000 psi and 8000 psi. The effect on the performance due to 
varying cowl lengths and flow Mach numbers was also investigated. The 
location of the impingement point of the cowl shock was determined for 
various cowl lengths and special attention was paid to the case of shock 
cancellation at the body surface. This work was done in collaboration 
with Dr. J.-L. Cambier whose two-dimens ional/axisymmetric nonequilibrium 
codes were used and during the course of this study a detailed grid 
sensitivity analysis of the computer codes was done. The adaptive gird 
code SAGE developed by Dr . Venkatapathy and Dr. Davies was also used in 
order to more accurately obtain the impingement points. (3) In support 
of the absorption experiment to be located at station N3 of the 16" 
shock tunnel nozzle, the NO (y) system [(0,0) and (0,1) bands] was 
studied in absorption with both quasi-one dimensional and two- 
dimensional axisymmetric reacting flow analyses. The required depth of 
the optical access, to minimize the effects of the hot boundary layer on 



the measurements, was estimated. The NO(y) band system was studied in 
very high-resolution with and without splitting of the 2 Z- 2 U 
transitions. The entire analysis was done with the assumption of 
chemical nonequilibrium and thermal equilibrium using a spectral code 
NEQAIR developed by Dr. Ellis Whiting. (4) A similar analysis was done 
at the exit of the nozzle in order to characterize the chemical state of 
the gas. (5) The spectra code was modified in order to do the 
spectroscopic study of the contaminants in the shock tunnel. The 
spectroscopic data for the species Fe, Cr, Ni and Cu are now 
incorporated into the code and the code has been handed over to the 
experimentalists who plan on using it to study the contamination levels 
in the shock tunnel. (6) The calibration of the Pinckney static pressure 
probes was initiated using a nonequilibrium parabolized Navier-Stokes 
code . 
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AXIAL VARIATION OF SURFACE PRESSURE 
WEDGE TEST - ARC 16” SHOCK TUNNEL 
Driver Pressure = 272.1 atm, Enthalpy = 10.5 MJ/kg 
0 = 11', V„ = 4107 m/s, T_ = 1475.4 K, p„ = 0.089 atm 
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2D COMBUSTOR MODEL - IMPINGEMENT STUDY 

IviACH CONTOURS AT THE LOWER CORNER 

P = 408.2 atm, H - 6.8 MJ/kg, M = 12 
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COMBUSTOR MODEL - IMPINGEMENT STUDY 

VIACH CONTOURS AT THE LOWER CORNER , 

P = 408.2 , H = 10.5 MJ/kg, M = 15.5 
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IMPINGEMENT STUDY - 2D COMBUSTOR MODEL 
Driver Pressure = 408.2 atm 
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GRID SENSITIVITY ANALYSIS - LOWER WALL PRESSURE 
IMPINGEMENT STUDY - 2D COMBUSTOR MODEL 
Mach 15.5 Case, Turbulent Flow, Driver Pressure = 408.2 atm 
p„ = 11.2 kPa, V„ = 4097.6 m/s, = 1412 K 
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GRID SENSITIVITY ANALYSIS - LOWER WALL PRESSURE (CORNER DETAIL) 
IMPINGEMENT STUDY - 2D COMBUSTOR MODEL 
Mach 15.5 Case, Turbulent Flow, Driver Pressure = 408.2 atm 
11.2 kPa, V„ = 4097.6 m/s, T„ = 1412 K 



Axial Distance - x(in) 









Pressure (atm) 




(ui) A - |jeM uioqoq aqi uiojj 3dub^stq 


Mach Number 



GRID SENSITIVITY ANALYSIS - TEMPERATURE PROFILE AT x = 40' 
IMPINGEMENT STUDY - 2D COMBUSTOR MODEL 
Mach 15.5 Case, Turbulent Flow, Driver Pressure = 408.2 atm 
p„ = 11.2 kPa, V„ = 4097.6 m/s, T„ = 1412 K 
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IMPINGEMENT STUDY - 2D COMBUSTOR MODEL 
MACH CONTOURS FOR TURBULENT FLOW 
M = 12, P = 544 atm, H = 6.6 MJ/kg 
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ARC 16" SHOCK TUNNEL - NOZZLE DIAGNOSTICS 
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During the last work period, I worked on the following projects: 

1 . I participated in the Ames Technical Paper Contest for Women. I presented work on hypersonic 
drag devices. I was one of the winners of this contest and will be attending the Society of Women 
Engineers Conference in June to present my winning paper. Also, a paper on this material has been 
submitted to the Journal of Spacecraft and Rockets for publication. 

2. I have also been working on computing the flowfield around a lunar return aerobrake concept 
proposed by Mike Tauber. I am working with E. Venketapathy on this project. Thus far an 
axisymmetric study of the aerobrake design has been completed for an aerobrake with flares at angles 
of 50°, 55°, and 60°. Three-dimensional calculations have also been completed at the same flare 
angles. If has been found that because of the large separated region created by the flares, the flow 
about the aerobrake is very complicated. Further investigation of this flow is anticipated and will 
include calculation of the base flow region. Comparisons between the axisymmetric solutions and the 
3-D solutions have been made for surface pressure, size of the separated region, and total drag 
produced. Validation of the drag calculations and general flowfield character are to be completed by 
calculating a similar flowfield that has been studied experimentally. A representative picture of this 
work is shown in Figure 1 . 

“3. I am continuing as point of contact for the University of Maryland. Maryland is involved in the CFD 
support effort for the Ames 16" shock tunnel. Thus far J.L. Cambier's axisymmetric TVD code 
(Emozart_2dA.f) has been transfered to the University of Maryland as well as the configuration of the 
ARC 16" shock tunnel nozzle and a generic combustor design. Cambier's 3-D code will eventually be 
transfered to the University of Maryland so that they can perform a full 3-D calculation of the test 
section flow with the combustor model in the test section. This study is being performed in order to 
determine the likelihood of choking the tunnel with the combustor model. 


4. My involvement with the CFD support of the 16" shock tunnel has increased in the last six months 



and will increase further when D. Prabhu leaves the project. I am also working with J.L. Cambier on this 
project A fully turbulent, nonequilibrium chemistry calculation of the nozzle flow has been completed. 
Information from this solution was used to characterize the flow in the test section and also to 
determine the optical character of the flow in the nozzle. A study of grid effects on the start-up 
solution has also been performed. This study was done assuming an inviscid, non-reacting flow. This 
start-up process begins at the breaking of the the diaphragm between the driven tube and the nozzle 
throat. A representative figure demonstrating grid effects is shown in Figure 2. A study of the entire 
start up procedure for a "best case” grid was also completed. This was also for inviscid, non-reacting 
flow. Figure 3. shows a time history of the start-up procedure. 

A study of the applicability of calculating quasi-one-dimensional (Q-1D) flow of a reacting gas in order 
to simulate the full axisymmetric, turbulent solution was performed. It is shown that by decreasing the 
nozzle radius by 13% (an 87% nozzle), the flow in the turbulent, reacting case could be simulated 
successfully at a fraction of the CPU time required for the axisymmetric calculation. The results of this 
study are shown in Figure 4. Note that the Q-1 D solution gives no information about boundary layer 
thickness or core flow profiles. This study will be continued further in order to determine if some "rules 
of thumb" exist for simulating axisymmetric, turbulent flows with finite rate chemistry with a Q-1 D code 
with finite rate chemistry. The Q-1D code is also being used to study the effects of varying the 
chemical reaction rates for the reactions taking place in the nozzle. 

Another study being done at this time is the simulation of the full driver tube, driven tube and nozzle 
flows in a time accurate manner. This simulation is axisymmetric, turbulent, and has finite rate 
chemistry. A representative figure of this flow is shown in Figure 5. In this figure, the diaphragm 
between the driver tube and the driven tube has partially ruptured and a shock has formed in the 
driven tube. 

m— «* > 

The results completed thus far for the Ames 16" shock tunnel will be presented at the AIAA 
Thermophysics Conference in July. 
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ARC 16" SHOCK TUNNEL - PROFILES AT NOZZLE EXIT 

A/A (geometric) = 169.8, R = 49.53 cm 
Shock Speed = 3.0 km/s. Reservoir Pressure = 420.5 atm 
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Figure 5 

Temperature Contours of Driver - Driven Section 
with Diaphragm Partially Open 
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ABSTRACT 

The single-expansion-ramp-nozzle (SERN) exper- 
iment underway at NASA Ames Research Center simu- 
lates the National Aerospace Plane propulsive jet-plume 
flow. Recently, limited experimental data has become 
available from an experiment with a generic noz- 
zle/afterbody model in a hypersonic wind tunnel. The 
present paper presents full, three-dimensional solutions 
obtained with the implicit Navier-Stokes solver, FL3D, 
for the baseline model and a version of the model with 
side extensions. Analysis of the computed flow clearly 
shows the complex 3-D nature of the flow, critical flow 
features, and the effect of side extensions on the plume 
flow development. Flow schematics appropriate for the 
conditions tested are presented for the baseline model and 
the model with side extensions. The computed results 
show excellent agreement with experimental shadow- 
graph and with surface pressure measurements. The 
computed and experimental surface oil-flows show the 
same features but may be improved by appropriate turbu- 
lence modeling. 


I. INTRODUCTION TO THE 
SERN EXPERIMENT 

Hypersonic air-breathing vehicles, such as the pro- 
posed National Aerospace Plane (NASP), require an 
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efficient integration of propulsion and aerodynamic sys- 
tems. A major element of the integrated design is the in- 
teraction between the jet-plume flow, the external flow, 
and the afterbody. The afterbody is designed as a single- 
(or one-sided) cxpansion-ramp-nozzle (SERN) to mini- 
mize friction drag while extracting thrust from the high 
pressure flow on the afterbody. This complex interaction 
can significantly alter the thrust, aerodynamic stability 
and overall performance of the NASP. 

Existing ground-based test facilities can only simu- 
late some aspects of the hypersonic flight environment. 
Therefore, accurate and validated computational fluid 
dynamics (CFD) codes are needed to provide simulations 
of realistic flight conditions. The experimental validation 
data should, of course, model as many of the realistic 
flight conditions as possible. The validated codes can then 
be used to provide the most reliable predictions of the in- 
crements in performance or design parameters associated 
with the differences between the available test conditions 
and the flight environment. 

An experiment has been conducted at NASA Ames 
Research Center which will allow greater insight into 
nozzle/afterbody flows relevant to the NASP. The SERN 
experiment was conducted in the NASA Ames 3.5 Ft. 
Hypersonic Wind Tunnel. This is a closed circuit, blow- 
down wind tunnel which has interchangeable, contoured, 
axisymmetric nozzles upstream of the test- section. 
Figures 1 and 2 illustrate the primary features of the 
SERN model. The side view of the model is a parallelo- 
gram. The leading- and trailing-edge angles of the model 
are 20°. The upper surface of the model forward of the 
cowl is a flat plate with a sharp leading-edge. A remov- 
able boundary-layer trip is located 4 in. downstream of 
tiie leading-edge. A perforated plate is located at the en- 
trance of the supply pipe to the plenum, followed by two 
screens. The internal-nozzle downstream of the plenum, 
was designed to provide uniform flow at the simulated 
combustor exit station. A nozzle cowl, 4 in. length, be- 
gins at the combustor exit station. The internal-nozzle 
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exit referred to in this study corresponds to the cowl exit 
station. 

Earlier versions of the model design included 15.2 
cm. wide extensions that would attach to the side of the 
model and effectively increase the ramp widlh. The en- 
tire model is immersed in a freestream flow of Moo = 7.3 
and high pressure air exits the nozzle with a jet to 
freestream static pressure ratio of approximately 11. 
Table 1 gives the flow conditions used for the present 
computations. A number of parametric studies were per- 
formed during the experimental investigation to under- 
stand the flow-field development under different operat- 
ing conditions. 

Figure 3 shows the basic flow features associated 
with the symmetry plane of the SERN experimental 
model, and the details are described below. The flow is 
three-dimensional, but most of the flow development can 
best be explained with a side view. Cross-flow views are 
given in a later section. The experimental model is ori- 
ented at a negative angle of attack with respect to the 
freestream. Leading-edge, or forebody, shocks are pre- 
sent both on the upper and lower surfaces of the model. 
Near the leading-edge, these shocks are limited to the up- 
per and the lower surfaces. However, away from the 
leading-edge, the flow on the sides of the model is influ- 
enced by these forebody shocks. The internal -nozzle flow 
development is dictated by the contour of this nozzle 
which is designed to be shock free. The internal-nozzle 
shape is flat on the upper and side walls and only the 
lower wall is contoured. Due to the one-sided, contoured 
wall design, the flow at the exit of the nozzle will vary in 
the vertical direction but will be two-dimensional. 

The interaction of the jet with the external flow is 
of main interest here. Similar to the plume flow associ- 
ated with rocket nozzles 1 * 2 the interaction between the 
expanding plume and the external flow produces an ouler 
plume shock, a shear-layer, and an inner plume shock 
(i.e., barrel shock). Since the ramp surface limits the 
plume flow development, the outer plume shock, the 
shear-layer and the inner plume shock will interact with 
the viscous wall layer on the ramp surface. The pressure 
differences between the upper and the lower surfaces of 
the model will induce a cross-flow from the lower to the 
upper side and a complex vortical structure will exist. 
The interaction between the plume flow and this side- 
edge flow will depend on the model geometry and flow 
conditions, such as angle of attack, etc. The static 
pressure ratio between the nozzle exit and the external 
flow primarily determines the size of the plume and the 
strength of the various flow structures associated vvilh the 
plume. The pressure and skin friction on the ramp 
surface determine the nozzle efficiency and the thrust. 
The pressure and skin friction values on the ramp surface 
are a result of the various interactions. In addition to the 
flow details mentioned above, the boundary-layers and 
shear-layers are expected to be turbulent. The flow may 
separate on the upper cowl surface due to adverse 
streamwise pressure gradient and may also be a plume 


induced cross-flow separation which will be discussed in 
the results sections 3 . The major plume features will be 
present whether or not the flow separates on the cowl 
upper surface. 

A variety of flow- field measurements have been 
taken both on the model surface and in the surrounding 
flow. At present, surface oil-flow patterns, shadowgraph 
visualization photographs and surface centerline static 
pressures are available and are presented in this paper. 
Complete surveys of the plume flow-field and skin fric- 
tion data will be presented in a subsequent publication 
which also describes, in detail, the experimental 
program. 

The experimental development has been supported 
by a computational effort. Earlier numerical studies 
were also used to help design tins experiment. These ef- 
forts are described in greater detail by Ruffin et. al. 4 and 
by Venkatapathy et. al. 5 However, these preliminary 
computations were performed on an earlier model design 
with side extensions. In the present paper, 3-D Navier- 
Stokes solutions and 2-D symmetry plane solutions are 
presented for the baseline model at the nominal test con- 
ditions. Critical flow features are identified and the re- 
sults are compared with experimental shadowgraphs, oil- 
flow visualization, and with surface pressure data along 
the model centerline. For the 3-D cases patched, and 
in some cases, solution adapted grids are used to map the 
complex 3-D geometry. Simulations of a version of the 
model with 15.2 cm. ramp side extensions are also pre- 
sented. The effects of side extensions and other 
geometric features on SERN flow structure are discussed. 

II. SOLUTION ALGORITHM 

The computations in this study are performed with 
FL3D, a three-dimensional, implicit time marching 
Navier-Stokes solver. This code has previously been 
applied to axisymmetric and three-dimensional generic 
rocket nozzle plume flows. L 2 All of the present calcula- 
tions is for laminar flow of a perfect gas. The perfect gas 
assumption is valid because both the freestream and jet 
static temperatures are sufficiently low. 

The numerical method is a LDU-ADI scheme with 
Roe’s averaging and MUSCL differencing. This numeri- 
cal formulation is detailed by Obayashi. 6 High order 
spatial accuracy is achieved by constructing MUSCL dif- 
ferencing with a differentiable limiter. The ADI sweeps 
of this formulation are quite efficient. The LDU-ADI al- 
gorithm is a diagonal algorithm requiring minimal CPU 
per iteration and is applicable to steady flows. 
Venkatapathy and Feiereisen 2 used tills method to predict 
plume flows accurately. 

This scheme has been applied to complex plume 
flows and the solutions were validated in terms of shock 
reflection location and plume structure for highly under- 
expanded axisymmetric plume flows. The calculations 
were performed both with supersonic external flow and 
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with quiescent ambient conditions. 4 The success of the 
solver in these limited validations permits confidence in 
the present solutions for the experimental SERN flow- 
field. 

Grid Generation 

The computational grids used in this study are cre- 
ated with an algebraic grid generation code 7 . This code 
utilizes bi-linear interpolation to provide smoothly vary- 
ing grid planes with nearly orthogonal grid cells. In 
order to map the complex SERN geometry, multiple grid 
zones are generated. The surface grid and the boundaries 
of each grid zone are shown in Figure 4 for the baseline 
model. This grid consists of 13 grid zones. Although not 
shown in Figure 4, each of the grid zones overlaps its 
neighbor by one common grid plane. The flow solver 
alternates grid zones on which it computes and applies the 
appropriate boundary condition on each of the intersect- 
ing grid planes. The computed results show seamless 
contours of all flow variables. The 3-D flow computation 
on the baseline model was performed with 338,000 total 
grid points and required approximately 5 hours of CPU 
time on a Cray YMP computer. 

Boundary Conditions 

The boundary conditions encountered in the pre- 
sent nozzle/afterbody calculations are inflow or outflow, 
wall boundaries, symmetry and reflection boundaries. 
All the boundary conditions were applied explicitly. At 
die supersonic inflow boundary, conservative flow vari- 
ables were frozen, and at the outflow, flow variables 
were extrapolated from the interior. At the wall bound- 
ary viscous, no-slip, adiabatic wall boundary conditions 
were applied. In the case of symmetry boundaries, the 
flow variables were simply reflected from the interior. 
The various boundary conditions were unified into one 
generalized routine so that any specified part of the 
computational domain could be updated with the proper 
boundary condition. This unified boundary condition 
formulation greatly facilitated the use of 3-D zonal grids 
needed to map the experimental geometry. 


III. 3-D RESULTS ON THE MODEL 
WITH SIDE EXTENSIONS 

One of the afterbody design considerations for hy- 
personic vehicles is the width of the afterbody surface. 
This geometric factor influences total vehicle thrust due 
to nozzle plume pressure and skin friction drag on the 
afterbody. To study the effects of ramp width, calcula- 
tions of a version of the model with 15.2 cm. side exten- 
sions were conducted. These calculations were per- 
formed at condition 1 shown on Table 1. For this case, 
the internal-nozzle flow was not computed. Instead, the 
jet flow at the nozzle exit is modeled by specifying invis- 
cid conditions and approximate boundary-layer profiles. 
The boundary-layer profiles are specified by using the 


Pohlhausen polynomial approximation 8 for the velocities 
and the Crocco-Busemann relation9 for compressible 
boundary-layers for the temperature profile. Figure 5 
shows the computational model for this case. For this 
simulation a solution adapted grid was used. The 3-D 
adapted grid package used is described by Davies and 
Venkatapathy 10 and the resulting adapted grid in the 
plume region is shown in Figure 6. The grid-adaption 
package clusters existing grid points in regions of high 
flow gradients. For the present case, the grid was adapted 
on gradients of both Mach number and static 
temperature. Figure 7 gives computed Mach contours in 
the symmetry and outflow planes. The solution was 
adapted on gradients of both temperature and Mach 
number so that good resolution of the plume shocks and 
shear-layers are achieved. 

The symmetry plane Mach contours in Figure 7 
predict the outer plume shock, shear-layer and inner 
plume shock (i.e., barrel shock) shown in the schematic in 
Figure 3. This calculation is performed for the model at 
zero angle of attack. In the outflow plane at the end of the 
ramp, the Mach contours also show these features and 
indicate that the width of the plume extends laterally to 
near the side-edge of the model. The Mach contours in 
this outflow plane also indicate a complex viscous struc- 
ture near the intersection of these plume shocks with the 
ramp and another complex structure just outside the 
plume boundary. These features can be further under- 
stood by studying Figures 8 and 9. Figure 8 shows simu- 
lated oil-flow and the Mach contours in the outflow plane, 
and Figure 9 shows particle traces and simulated surface 
oil-flow pattern. For the model with side extensions, the 
freestream flow expands in the region above the side ex- 
tension ramp surface. The magenta particle traces in 
Figure 9 correspond to a large vortex which is generated 
in this expanding region above the side extensions. As the 
flow proceeds downstream, this side extension vortex 
lifts off the ramp surface, and follows the outside edge of 
the outer plume shock. The blue particle traces 
correspond to a vortex which begins from the trailing- 
edge comer of the cowl and follows the outer boundary 
of the plume along the surface of the ramp. The most 
prominent features in the predicted oil-flow pattern are 
cross-flow separation lines which extend from the cowl to 
the ramp trailing-edge. These separation lines are clearly 
within the plume external shock boundary and approach 
the symmetry line of the model as they proceed 
downstream. These features correspond to cross-flow 
separation induced by the inner plume shock (barrel 
shock) within the plume. The expanding cross-flow from 
the symmetry line experiences a strongly adverse 
pressure gradient from the internal plume shock and 
cross-flow separation is induced. 

Based on the results presented, a schematic of the 
cross-sectional SERN flow is developed and is shown in 
Figure 18b. This flow schematic is appropriate in the 
plume region of the flow at the conditions given in Table 
1. The large side extension vortex lifts off tire ramp sur- 
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face and interacts with the outer plume shock. Near the 
ramp, the inner plume shock induces a primary cross- 
flow separation vortex whose separation line is very 
prominent in the surface oil-flow pattern. Underneath 
this vortex, a secondary cross-flow vortex is found to ex- 
ist. 


IV. RESULTS ON THE BASELINE MODEL 

Computations on the baseline SERN model without 
side extensions were conducted at condition 2 shown on 
Table 1. For this simulation the internal -nozzle flow was 
computed and solved coupled to the other grid regions. 
Predicted Mach contours in the symmetry and outflow 
planes for the baseline experimental model are shown in 
Figure 10. The symmetry plane shows the same features 
in the plume region as previous computations on the 
model with side extensions. The calculations on the base- 
line model and the experiment were performed with the 
model at a -1.1° angle relative to the freestream and the 
forebody produces a shock on the upper side of the 
model. From the outflow plane we see that the width of 
the plume extends to the side-edge of the body. Particle 
traces and velocity vectors which are not shown indicate a 
number of vortices near the model side-edge. The 
overall plume flow structure in the cross-flow plane for 
the baseline model is shown in Figure 1 8a for these flow 
conditions. The outer plume shock extends laterally over 
the side-edge of the model and approaches the side comer 
of the model near the ramp. This shock induces a vortex 
on the side of the model. As in the model with side exten- 
sions, the inner plume shock induces a primary cross- 
flow separation vortex, which is accompanied by a 
smaller secondary vortex. Because the baseline model 
has no side extensions, the large side extension vortex 
present in the previous results does not exist for the 
baseline model. 

Surface pressure measurements have been taken at 
a variety of locations on the model ramp. At present only 
data along the centerline is available. Figure 1 1 shows a 
comparison of the experimental and computed static pres- 
sure along the ramp centerline. The computed results 
show excellent agreement with the experimental data. 

The symmetry plane flow is of interest for two rea- 
sons. First, although experimental surveys of several 
cross-flow planes will be available, the most detailed data 
will be obtained along the symmetry plane. Also, very 
fine grids can be used in 2-D simulations to resolve criti- 
cal flow features without prohibitive CPU requirements. 
For these two reasons, it is of interest to know if the 
three-dimensional features affect the symmetry plane 
flow. To study 3-D effects on the symmetry plane flow, 
2-D symmetry plane results were compared to the pre- 
dicted symmetry plane flow from a 3-D calculation. The 
symmetry plane grid contains 80 x 60 grid points. Figure 
12 shows Mach number surveys at three stream wise loca- 
tions for the 3-D and 2-D cases on the same grid. The 
x=24.8 in. station is upstream of the nozzle cowl. The 


x=41.9 in. station is in the plume region midway between 
the internal-nozzle exit and the end of the ramp, and the 
x=51.1 in station is in the plume region at the end of the 
afterbody ramp. It is found that the 3-D features do not 
significantly affect the symmetry plane flow until near 
the trailing-edge of the model. Similar comparisons of 2- 
D and 3-D results in the symmetry plane 4 also revealed 
little difference between the two simulations. However, 
overall the flow does expand somewhat faster in the 3-D 
case due to pressure relief in the cross-flow direction. 

Calibration data obtained with the Mach 7.3 wind- 
tunnel nozzle used during the experimental investigation 
show the presence of a moderate degree of expansion 
(source flow) in the test-section. For most of the compu- 
tational simulations, this gradient was neglected, and uni- 
form external-flow conditions were imposed on the in- 
flow boundaries of the computational domain. To evalu- 
ate the effect of the non-uniform external flow and to 
study any interaction between the non-uniform free- 
stream and the jet plume, numerical experiments were 
conducted with uniform and expanding free-stream flow 
conditions. These numerical experiments were limited to 
2-D, symmetry plane computations with a 159*1 19 grid. 
In the tunnel test-section, a normalized pressure change 
of Ap/p^ = -0.138 per meter and Mach number change of 

4M m = 0.178 per meter were measured. In the 2-D 

simulations, these gradients were imposed and through 
isentropic relationships! 1 the freestream conditions and 
the conservative flow quantities at all the inflow and top 
boundary grid points were computed. In the uniform 
freestream case, the flow conditions corresponding to 
condition 2 in Table 1 was maintained at all grid points on 
the top and inflow boundary surfaces. Both cases were 
computed with the same nozzle sonic throat conditions 
for the internal-nozzle. Pressure contours from the com- 
puted solutions with uniform and non-uniform 
freestream conditions are shown in Figure 13 (a) and (b) 
respectively. Static pressure profiles, along the z-direc- 
tion, at three streamwise locations, are compared in 
Figure 14 (a), (b) and (c). The computed solutions show 
that the non-uniform freestream conditions affect only 
the flow outside the jet-plume region. The upper fore- 
body shock shape and the conditions behind the shock are 
affected by the conditions upstream of this shock, but the 
shock acts as a buffer to reduce the effect of the non-uni- 
formities in the free-stream. The outer plume shock has a 
similar damping effect, and the flow-field in the vicinity 
of the jet closely resembles that of the uniform flow case. 
The ramp surface pressure and the flow inside the plume 
arc unaffected by the non-uniform external flow. 
Though the numerical experiments were carried out only 
along the symmetry plane, as a 2-D study, these conclu- 
sions are also expected to hold in the 3-D case. 

Symmetry plane Mach contours on a 2-D 300x300 
fine grid are shown in Figure 15. This fine grid simula- 
tion resolves the outer plume shock, shear-layer, and in- 
ner plume shock quite well. It also indicates that the 
downward curve of the upper cowl surface induces a 
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edge of the cowl. The curved outer boundary of the cowl 
followed by the shear-layer between the under-expanded 
jet and the freestream produces an effective compression 
corner. 

An experimental shadowgraph photograph and a 
simulated shadowgraph based on the fine grid symmetry 
plane calculation are shown on Figure 16. Again wc sec 
the forebody shock, boundary-layer separation near the 
cowl trailing-edge, outer plume shock, shear-layer, and 
barrel-shock. The location of these features agrees very 
well with the computed results. One difference, how- 
ever, is that the experimental boundary -layers and shear- 
layers arc turbulent. The experimental shadowgraph 
shows a thick turbulent boundary-layer on the model 
forebody which feeds into a thick plume shear-layer. The 
laminar flow calculations of course do not reproduce 
these turbulent details but predict the overall symmetry 
plane flow quite well. 

Predicted surface oil-flow patterns on the after- 
body ramp for the baseline model and the model with side 
extensions are shown in Figures 17b and 17c. The exper- 
imental pattern is shown in Figure 17a. In the 
experimental oil pattern, the flow appears to be two- 
dimensional inside the cowl, with no separation at the 
corner. The pattern on the ramp is symmetrical and the 
centerline flow is straight down the ramp. The surface 
flow direction turns increasingly outboard with 
increasing distance from the centerline. The prominent 
cross-flow separation lines on the ramp predicted in both 
of the computed results arc also observed experimentally. 
As is evident from the oil-flow patterns and indicated in 
schematics in Figures 18a and b, the primary cross-flow 
separation lines are found to exist closer to the symmetry 
line for the baseline model relative to Lhe model with side 
extensions. The rale at which these separation lines 
approach the symmetry line toward the trailing-edge is 
different in the predicted vs the experimental results. 
The details of this viscous dominated effect may be 
strongly influenced by turbulence which is not modeled, 
in all of the oil-flow patterns a secondary cross- flow 
separation line is also seen near the side-edge of the 
model. In addition, the experimental oil -flow pattern 
shows a streamwise separation line on the upper surface 
of the cowl near the trailing-edge. This separation is also 
apparent from the computational and experimental 
shadowgraphs in Figure 16. 


CONCLUDING REMARKS 
The latest version of FL3D, an implicit, 3-D 
Navier-Slokes solver is used to predict the experimental 
hypersonic flow-field. Multiple, and in some cases, solu- 
tion adapted grids arc used to map the 3-D noz- 
zle/afterbody geometry. TVD enhancements are used to 
model critical flow features. These modifications have 
greatly improved the efficiency of the code and its ability 
to provide accurate modeling of complex flow-fields. 
Three-dimensional calculations of the two SERN geome- 


tries and 2-D symmetry plane calculations have been per- 
formed which simulate the cold-jet perfect gas experi- 
mental conditions. Flow schematics appropriate for die 
conditions tested arc presented for the baseline model and 
the model with side extensions. 

The computed results show excellent agreement 
with experimental shadowgraph and with centerline sur- 
face pressure measurements. The computed and experi- 
mental surface oil-flows show favorable agreement but 
demonstrate the need for turbulence modeling. At pre- 
sent, several candidate algebraic and two equation turbu- 
lence models arc being validated in jet plume flows. 
Turbulent code validation cases will be conducted and 
presented along with detailed calibrated experimental 
data in a subsequent publication. All of the How features 
observed experimentally arc also seen in the computed 
results, thus the success of the current predictions with 
the available data is quite promising. 


Table 1 . 



Condition- 1 

Condition-2 

Moo 

7.3 

7.4 

a 

0.0 

-1.1 

Too (°K) 

57.2 

71.5 

Poo 

756 Pa 

673 Pa 

Rc/m 

1.2x 10 8 

1.5x 10 s 

Ptjcl/Poo 

300 

337 

TtjeAco 

5.0 

3.9 

PjCl/Pco 

10.6 

= 10.7 

Mjct 

2.83 

- 2.87 

Note: Both 

the freestream and jet gas 

arc air. 
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Figure 1. Schematic of the baseline SERN Model 



Figure 2. Photograph of the baseline SERN Model 



M = 7.3 



Figure 3. Symmetry plane schematic of the SERN flow-field. 
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Figure 4. Surface grid and grid zone boundaries for baseline model computations. 



Figure 5. Computational body for model with 15.2 cm. wide side extensions. 
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Figure 6. Solution adapted grid for computations on model 
with side extensions. 



Figure 8. Simulated surface oil flow pattern and Mach 
contours in uie ouuiow plane lor model witn side extensions. 



Figure 7. Computed Mach contours in the symmetry and 
outflow planes in the plume region or the model with side 
extensions. Solution adapted grid is used 



Figure 9. Simulated surface oil flow pattern and particle 
traces lor the model with side extensions. 



Figure 10. Computed Mach con tours in the symmetry and 
outflow’ planes for baseline model. 
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Figure 11. Comparison of computed and experimental surface static pressure along the baseline model centerline. 



Figure 12. Comparison of symmetry plane Mach number at three streamwise stations for 3-D and 2-D calcula- 
tions. Solid line = 3-D solution, Dashed line = 2-D solution on the same symmetry plane grid. 





Figure 13. Comparison of symmetry plane pressure at three streamwise stations for flow with uniform freestream 
and flow with expanding freestream variation. Solid line = uniform freestream, Dashed line = expanding 
freestream. 













Figure 14. Comparison of symmetry plane pressure contours for a) flow with uniform freestream 
and b) with expanding freestream variations 



Figure 15 Computed Mach contours in the sy mmetry plane of the base line model 
using a 2*D fine grid. 



Figure 16. a) Experimental and b) computational shadowgraphs in the symmetry plane of the 
baseline model Simulated shadowgraph is based on 2-D fine grid results. 
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Figure 17. Experimental and simulated oil-flow on the ramp surface. 
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Figure 18. Cross-sectional schematics of SERN flow field in the plume region. 
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Abstract 

The axis singularity problem for three dimen- 
sional blunt body flows is addressed in this pa- 
per. Two methods of eliminating the axis-singularity 
difficulty are outlined for use with finite-difference 
codes. In the first method, the three-dimensional 
Navier-Stokes equations are reformulated using a re- 
defined grid Jacobian that is non-zero at the sin- 
gular line. The resulting equations along with ap- 
propriate boundary procedure at the axis are solved 
using the finite-difference form of a flux-split, up- 
wind scheme. In the second method, the three- 
dimensional Navier-Stokes equations are solved us- 
ing the Roe flux-difference splitting method with 
a finite- volume based method for evaluation of the 
metrics and grid Jacobian, and appropriate bound- 
ary conditions. Solution symmetry is preserved in- 
dependent of the orientation and location of the sin- 
gular line boundary with respect to the free-stream 
velocity vector. Inviscid solutions are compared to 
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the results of Lyubimov and Rusanov 1 who em- 
ployed a shock-fitting technique. Viscous and real- 
gas computations are also presented demonstrating 
the generality of the two techniques. 

Introduction 

Most three dimensional (3-D) grids constructed 
to solve for flow over a blunt body consist of series of 
grid planes distributed radially about a common or 
singular line. Generally, this singular line extends 
from the nose of the body. It is referred to as a 
singular line because it introduces a mathematical 
singularity into the solution algorithm, which can 
cause non-physical perturbations in the flow solu- 
tion. Dealing with the singular line is one of the 
most vexing problems in the computation of 3-D 
blunt body flows. A number of recent references 2-6 
have reported the occurrence of non-physical solu- 
tions in the vicinity of the singular line. 

Generally some form of extrapolation is used 
to determine the Jacobian, metric terms, and flow 
quantities along the singular line. The most com- 
monly used methods include a simple zeroth-or first- 
order extrapolation-averaging boundary condition. 
These methods can cause numerical problems in the 
flow solution. Irregularities are created in the pres- 
sure and velocity field.This is readily apparent if a 
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grid where the singular line coincides with the stag- 
nation line of the flow is used. A zeroth-order ex- 
trapolation causes the peak pressure to occur away 
from the stagnation point, and the radial velocity 
along the stagnation line will not be zero. Simi- 
lar behavior has been observed both by Peery and 
Imlay 2 , and Josyula and Shang 4 . 

This study attacks the singular line boundary 
condition problem with two different approaches. 
The first is a reformulation of the governing equa- 
tions based on a cylindrical-type grid Jacobian that 
is non-zero on the singular line. This allows the 
flux emanating from the singular line to be prop- 
erly computed. This method in essence moves the 
singularity to two source-like vectors that are not 
evaluated along the singular line. In the second ap- 
proach the governing equations in generalized coor- 
dinate system, obtained from the the Cartesian sys- 
tem, are solved using Roe's approximate Riemann 
solver along with MUSCL extrapolation scheme. A 
finite-volume based formulation is utilized in evalu- 
ating the grid Jacobian and the metric terms. With 
second order accurate boundary procedures, the sin- 
gular line problem is eliminated in this approach. 

Because the singular line problem is due to 
fluid dynamics and the evaluation of metric and flux 
terms, the derivation of the techniques and most of 
the solutions shown in this work are for perfect gas 
^ etwo approaches are tested by computing a 
Mach 20 mviscid flow over a sphere for comparison 
with accurate solutions reported in Ref.l. To eval- 
uate the applicability of the solution techniques for 
a wide range of problems, a grid refinement study 
along with numerous computations for viscous, adi- 
abatic and isothermal walls, three-dimensional ideal 
gas conditions and three- dimensional reacting-gas 
conditions were conducted. 


Method One: Reformulation 
of the Governing Equations 

A typical grid over a three dimensional blunt 
body is shown in Fig. 1. It consists of a series of 
radial grid planes that meet along a common or sin- 
gular line usually emanating from the nose of the 
body. Cell volumes adjacent to the singular line are 
wedge-shaped instead of cube-shaped. In the follow- 


ing discussions the singular line is assumed to lie on 
the x-axis. 

The solution procedure normally employed by 
CFD algorithms is to transform the governing equa- 
tions and grid coordinates from a Cartesian (x,y,z) 
coordinate system to a generalized (£, 77, Q coordi- 
nate system. In this and subseqent discussions in 
this section f-lines march along the longitudinal- 
, »?-lines the circumferential-, and C-lines the body 
normal-direction. The Navier-Stokes equations ex- 
pressed in vector form in generalized coordinates are 


dQ dE dF dG dR OS dT 

dt da + a v + ac ~ at + at + dc (1) 

with 

Q = J-'Q (2) 

£ = ■ r ‘[f r £ + { v .F + {«cj 

f = J ~' [’)*£ + <hF + t),a] 

S = ■'-‘[g-e+c./’+gc] 

where E,F, and G are the inviscid Cartesian 
flux vectors. The grid Jacobian, J, used in the 
above expressions and to compute the metric terms, 
is defined as 


J 1 = Mvi* < - y< Z„) - x v (y ( z c - y (Z( ) 

+*<:(!/« - y„*e)] (3) 

Most modern CFD algorithms use upwind dif- 
ferencing, splitting the inviscid flux vectors into pos- 
itive and negative components according to the di- 
rection of signal propagation. The term ,f or in- 

*?“?- , is s ? Ut <"to f components. In 

tfle finite- difference formulation, the evaluation of 

W at the points adjacent to the singular line 
involves evaluating the E + flux at the singular line 
Along the singular line, the spatial derivatives 
in the circumferential direction, and are 

all zero making the grid Jacobian defined by eq. (3) 
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zero. This causes the E+ flux to be incorrectly set 
to zero along the singular line. Essentially, no in- 
formation is allowed to pass from the singular line 
boundary to the interior causing spurious perturba- 
tions in the computed flow solution. 

A finite value for the singular line Jacobian 
can be obtained by extrapolation from interior grid 
points or by setting y v and z v , or x ^ and y ^ depend- 
ing on the axis orientation, to small finite values. As 
will be shown these approximate methods yield un- 
satisfactory results and lead to nonphysical flowfield 
perturbations in the region of the singular line. 

The first method for treating the singular line 
boundary problem presented in this study involves 
re-defining the grid Jacobian such that it is no longer 
zero on the singular line. The grid Jacobian is writ- 
ten as the product of two quantities 

J~ l = rJ~ l (4) 

where r is the radial distance from the singu- 
lar line, (y 2 + z 2 )i. The modified Jacobian, J _1 , 

is non-zero throughout the entire computational do- 
main. It is equivalent to the grid Jacobian obtained 
if the Navier- Stokes equations were originally ex- 
pressed in cylindrical, (x,r,0), coordinates and can 
be computed from the expression 

^ _1 = [ x d r v G c - T a 8v) - x v( r dc - T a e i) 

+ x d r dr, - (5) 

If Eq. (4) is inserted in Eq. (1) and the chain 
rule used to remove the r-term from the partial dif- 
ferential equations, the resulting equation is: 

dQdEdFdGS 8R 8S df W , x 
dt + dt + dT) + dC + r - dt + dT) + d( + r (6) 

The vector quantities are now scaled by the 
modified, or cylindrical, grid Jacobian. 

Q = J~ l Q 

E = + + 

F = 1 rj x E + Tj r yF + tj 2 G 


G = J 1 Q x E + C y F + ( Z G 

The viscous vectors are similarly scaled. The 
metric quantities, and Cartesian 

flux vectors are unchanged. Equation (6) is identi- 
cal to that obtained if the cylindrical Navier-Stokes 
equations are transformed to generalized coordi- 
nates. The additional source-like vectors, H and W , 
are defined as 

H = [r { E + r v F + r c G] (7) 

W = [r ( j? + r v S + r c T ] 

The singularity has, in effect, been transfered to 
these two source-like vectors. Since the source vec- 
tors are not evaluated along the singular line where 
r = 0, the source vectors are never singular. For 
grids over bodies of rotation where x v = r v = 0, the 
source vector, H , is given by 

pv r 

puv r 

H = J -1 pvv r + p cosO (8) 

pwv r + p sinO 
. ( e + p)v r . 

v T — v cos9 + w sin6 
0 = tan -1 (— ) 

V 

For purposes of comparison, the same zeroth- 
order extrapolation/averaging boundary condition 
along the singular line was used with both the for- 
mulation using the customary Cartesian grid Jaco- 
bian and the new formulation using the cylindrical- 
type grid Jacobian. For zero angle of attack, the 
velocity components v and w are set to zero along 
the singular line which is also the stagnation line for 
these cases. For non-zero angle of attack cases, the 
Cartesian v velocity component is still zero due to 
symmetry. The w component is obtained by extrap- 
olation /averaging. 

Method Two: Three-Dimensional 
Roe Upwind Solver 

In the second method, both the Navier-Stokes 
equations and the Euler equations are solved using 
the upwind flux-difference splitting scheme due to 
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^"° e • The equations in generalized coordinates are 
obtained from the Navier-Stokes or Euler equations 
written in the Cartesian coordinate system. The 
transformation from the physical (Cartesian coor- 
dinate) to the computational (generalized) coordi- 
nate system is performed numerically. The trans- 
formed equations in the computational coordinate 
are then solved using Roe’s approximate Riemann 
solver. The singular line problem is easier to study 
with an explicit formulation rather than an implicit 
formulation. Thus, only an explicit upwind formu- 
lation is studied and reported in this paper. 

In the finite-difference representation, appropri- 
ate flux evaluation in the neighborhood of the sin- 
gular line becomes an issue since the transformation 
from the physical to the computational coordinate 
is singular at the singular line. Instead if a finite- 
volume representation is followed, then the metric 
and Jacobian are determinate and the singular line 
difficulty can be avoided. In the present method, 
following Vinoukur 9 and Obayashi 10 , we combine 
the finite-difference method of solving the governing 
equations with the finite- volume method of evaluat- 
ing the grid metrics and the Jacobian. For all the 
flows solved using the present method, a finite- dif- 
ference grid was selected such that grid points were 
located along the singular line (see Figure 1). 

The finite-difference representation of equation 
(1), with the upwind-biased flux-difference splitting 
of m viscid terms due to Roe 7 and with central- 
difference representation of the viscous terms is writ- 
ten as : 


|Q n+1 _ 0»1 _ 

A < [(^(i+£,i,A0 - E {i _^ jk) ) + (F(ij + I <k ) - F (ij _^ k) ) 

+ (G(i,j+i,k) ~ £(.,;-*,*)) - ~ 

- (T(ij <k+ l) - 7(,j,/fc_$))](9) 

In the above equation, the change in Q at the 
grid point (i,j,k) between two time levels is repre- 
sented as the the difference in the fluxes (inviscid and 
viscous) around the surfaces elements of the control 
volume placed at the grid point (i, j, k). These fluxes 
are evaluated at the cell interfaces placed half-way 
between the grid points in each computational di- 
rection. The inviscid fluxes are upwind differenced 


and are evaluated using Roe’s flux-difference split- 
ting method. For example, the inviscid flux in the 
^-direction at (*’+ 1/2, j, k) interface, utilizing Roe’s 
approximate Riemann solver, is given by 

E(i+i/2,j,k) = 1/2 [E(Qr, Sf +1/2iJi *)+ 

E(Ql, Si+l/2,j,k)~\MQL , Qr, S,- + i /2,j, k )\(Q r-Q l)} 

(i°) 

where |A| is the Roe-averaged Jacobian of E and Q L 
and Q R are the state variables to the left and right 
of the interface at (» + 1/2 J, k ) and Sf +1/2 is the 
surface area vector of the interface. The surface area 
vectors are related to the metric terms. For exam- 

p * e § L,k = J~ X j£r l + ZyJ + £ z kJ . Roe-averaged 

flux-difference representations, similar to equation 
(10) are used for inviscid fluxes in tj and ( direction. 
The viscous fluxes are central differenced and each 
viscous flux at the half- cell interfaces involves the 
state vector Q and the metric terms from the neigh- 
bouring grid points. For example, the viscous flux 
in the £ direction at the cell interface (i +1/2 i k) 
will be of the form ’ 


Ri+i/2,j,k = R{Qi,j,k,Qi, j+1 , k ,Q iJ _ 1 ' k , 

Q*,j,k+ 1 , Qi+l,j+l' k , Qi+l,j-i' k , 

Q*+i,j,k+uQi+ij,k-i , Sf +1/2> . >fc , s* +1/2jk , 


®t+l/2,i,*) ( 1J ) 

The above difference scheme can be spatially 
first, second or third order accurate depending on 

wL a £? rOXimati ° n Used for Ql ™* Qr- Through 
MU SCL extrapolation 7 for Q L and Q R in equation 
(10), second and third order spatial accuracy can be 
achieved. When applying MUSCL extrpolation, the 
primitive variables are limited so as to be Total Vari- 
ation Diminishing (TVD) and in the present study 
a differential limiter 7 is used. In addition, in evalu- 
ating the flux Jacobian matrix, |A|, entropy correc- 
tion is required when the magnitude of the eigenval- 
ues become small. The entropy correction used here 
is similar to the one suggested by Yee 11 and used 
y ayashi but with the coefficient term propor- 
tional to the second derivative of the pressure. The 
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entropy corrected eigenvalue A of the left and right 
states are 

A = (A 2 + c 2 )/2c 

*/ l-M < f (12) 

where 


€ “ C l(l^ Pl/|Pavcra5c|)(|^max|)* (13) 

In equation 13, p is the pressure and 0.2 < c\ < 1.0. 

The evaluation of inviscid and viscous fluxes in- 
volve not only the state vectors, Q, but also the sur- 
face area vectors ( the metric terms ) and the Ja- 
cobian. JThe^metric terms, represented by the area 
vectors , S v and S c and the Jacobian of the trans- 
formation (which is related to the volume by J~ l 
= where V is the volume) are all evaluated by 
constructing a finite- volume around each grid point 
(i j,k). For the grids used in the present study, the 
volumes are made up of 6-sided elements away from 
the singular line and 5-sided elements near the sin- 
gular line. At all interior points, where the volume 
is 6- sided , the surface area vectors and volume can 
be evaluated 9 as : 


Sij,k - g( ? (<,j+i,fc+i) - r(«,i-i,fc-i))x 


('(•j-i.fc+i) - *(« f i+i,*-i)) (14) 
where r is the position vector, and 
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V ” + + S ?,i+ i t * + S lj-i,k+ 


Sfj.fc+i + Sf iiifc _ 1 )«(r( i+liJ+1>fc+1 ) - ?(«-ij-i,fc_i))(15) 
The surface area vectors required in expressions 
similar to equations (10,11) are at the cell interfaces 
and are computed by averaging the area vectors at 
the grid points. 

^<+1/2, j,k = 2 (^+ 1 , j.fc + S lj,k) (l 6 ) 


computed by applying appropriate boundary condi- 
tions, evaluation of the metric terms and the Jaco- 
bian at the singular line are not required. But the 
solution at grid points adjacent to the singular line 
requires the surface area vectors half-way between 
these points and the singular line. The required sur- 
face area vectors and the volume elements, for all the 
grid points adjacent to the singular line, are evalu- 
ated by constructing appropriate surface and volume 
elements around these grid points instead of using 
equations 14, 15 and 16. These evaluations are very 
simple and are not shown. Once the appropriate 
area vectors and the Jacobian are evaluated at all 
the interior points, the state vector can be updated 
from one time level to the next. 

Using equation (9) the solution vector at all the 
interior points are explicitly updated. To compute 
axisymmetric flows, only three computational planes 
along the tj or circumferential direction are used and 
the solution is sought only on the mid-plane (see Fig- 
ure 9). The state vector on the adjacent planes are 
obtained through rotation. Three planes along the 
T] direction are required to impose axisymmtery. At 
the boundaries, such as wall, inflow, outflow or sym- 
metry , appropriate boundary condition procedures 
are employed to extrapolate the state vector Q to the 
boundary surfaces. Second order accurate extrapo- 
lation , involving Q(2J,k), Q(3 j,k) and Q(4 j,k) grid 
points result in Q(1 j,k) at the grid points on the sin- 
gular line boundary. 

For three-dimensional flows, a complete three- 
dimensional grid is used. If the flow is bi-symmteric 
then the grid is constructed in only one-half of the 
three-dimensional volume. The boundary conditions 
for the wall, inflow or outflow boundaries are similar 
to the axisymmteric (3-plane) boundary conditions. 
But for the grid points on the singular line, in addi- 
tion to the second-order accurate extrapolation, an 
averaging procedure (in the tj direction) at the sin- 
gular line for (p,p,\U\) is applied and the velocity 
|f/| is rotated to coincide with the velocity direction 
of the state vector at the adjacent grid point, (2 j,k). 
The accuracy of the present method is demonstrated 
through numerical examples in the next section. 


The flux-difference terms in the right hand side 
of equation 9 can be evaluated accurately without 
any difficulty at all interior points. Since the solution 
at the singular line boundary is not solved for but 


Results 
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Solutions computed using both methods are 
compared against data produced by Lyubimov and 



Rusanov 1 . They utilized a shock-fitting iterative 
difference scheme, referred to as a constant direc- 
tion scheme, to compute inviscid flow over simple 
bodies of revolution at various Mach numbers. The 
method used the non-conservative form of the gov- 
erning equations to solve for the vector of primitive 
variables X = [/>, u,v, u;,p] as well as the function F 
which defined the shape and location of the shock 
subject to the initial and boundary conditions. The 
results given in Ref. 1 are particularly useful for this 
study because they are presented as a series of ta- 
bles showing shock standoff distance and numerical 
values of normalized density, pressure, and velocity 
between the body and shock at various angles along 
the body surface. This allows a one-to-one com- 
parison between an independent calculation and the 
results computed in this study. 

The general test case selected was that of Mach 
20, inviscid perfect gas flow over a sphere at various 
angles of attack. This case was chosen to permit 
comparison with Ref. 1. The techniques presented 
are equally effective for viscous and real gas flows as 
will be demonstrated. 

Results Using the Cylindrical Formulation 

The first series of line plots concerns solutions 
for zero angle of attack, effectively axisymmetric flow 
over a sphere. Results are shown for both the 3- 
D Cartesian grid Jacobian and 3-D cylindrical grid 
Jacobian versions of the codes. Henceforth, the 
Cartesian grid Jacobian and cylindrical grid Jaco- 
bian solutions will be referred to as the Cartesian 
and cylindrical solution respectively. Both codes uti- 
lize shock-capturing with first-order accurate differ- 
encing used for the inviscid fluxes. Both codes used 
the same grid and the same zeroth-order boundary 
condition along the singular line. To give the grid Ja- 
cobian some non-zero but infinitesimally small value 
for the Cartesian grid Jacobian case the points on 
the singular line were given finite values of y and z 
corresponding to the points being distributed on a 
circle of radius 1.0 x 10 -7 normalized with respect 
to the sphere radius centered on the x-axis . The 
metric terms are identical for both cases. Thus, the 
only difference between the two codes is the use of 
the Cartesian or cylindrical grid Jacobian formula- 
tion. Data from Ref. 1 is also included on the charts 
for comparison as well as data from an axisymmetric 


version of the code. 

Figures 2a and 2b show normalized values of 
density and pressure along the body surface. The 
stagnation point and singular line correspond to 
0 = 0. The 3-D cylindrical, axisymmetric, and 
Ref. 1 data all compare very closely. The Carte- 
sian code, however, has difficulty near the singular 
line. The peak values of density and pressure occur 
away from the stagnation line and are too high. The 
radial velocity is not zero along the stagnation line 
as it should be. Away from the singular line, the re- 
sults from the Cartesian code gradually approaches 
that of the other solutions. A close-up of the surface 
pressure around the singular line region is shown in 
Fig. 2c. The Cartesian solution is clearly unsatis- 
factory exhibiting a sawtooth pattern. This is the 
consequence of the positive longitudinal flux being 
incorrectly set to zero along the singular line. The 
cylindrical code does a much better job of represent- 
ing that flux. Note also the effect of the zeroth-order 
boundary condition on the 3-D cylindrical solution. 
The axisymmetric solution does not have a grid line 
at y = 0 but instead has two symmetric grid lines 
above and below the x-axis. 

Figure 2d shows computed shock standoff dis- 
tance as a function of 6 around the body. Once 
again the cylindrical, axisymmetric, and Ref. 1 data 
correspond very closely. Another consequence of us- 
ing the Cartesian code is that the stagnation line 
shock standoff distance is underpredicted by more 
than 10%. Away from the singular line, the shock 
standoff distance is overpredicted. 

Normalized density, pressure, and u-velocity 
profiles along the stagnation line are shown in Figs. 
3a-c. These plots also highlight the fact the 3-D 
Cartesian code underpredicts the shock standoff dis- 
tance while the other three solutions are essentially 
identical. Figures 4 and 5 show density and pres- 
sure profiles along the 45 and 75 degree lines. The 
peak values computed by the two shock capturing 
codes are somewhat lower than the shock-fitted val- 
ues from Ref. 1 as would be expected. The Carte- 
sian code does somewhat better along these lines 
that are removed from the singular region. Pres- 
sure and shock standoff distance are somewhat over 
predicted. 

The cylindrical formulation is equally effective 
when applied to viscous flows. Mach 31.7 perfect 
gas flow was computed over the Aeroassist Flight 
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Experiment (AFE) vehicle, shown in Fig. 6, using 
both the Cartesian and cylindrical formulations with 
the same singular line boundary conditions used in 
the sphere calculations. The Cartesian code had a 
great deal of difficulty resolving velocity and pres- 
sure near the singular line. A surface pressure profile 
plot along the pitch plane is shown in Figure 7 . With 
the Cartesian solution, the distortion around the sin- 
gular line causes two localized pressure peaks away 
from the stagnation point. Clearly, this is a non- 
physical distortion caused by the numerics of the 
solution process. This pressure distortion affects the 
pressure distribution in the entire subsonic region of 
the flow, which for the AFE vehicle is most of the 
forebody. This illustrates how important it is to ac- 
curately resolve the flow near the singular line. The 
cylindrical solution shows a smooth pressure distri- 
bution with the peak value at the stagnation line. 

Finally, to show how readily extendable the 
cylindrical formulation is to real-gas flows, velocity 
magnitude contours using both Cartesian and cylin- 
drical formulations are presented in Figs. 8a and 8b 
for Mach 31.7 viscous flow over the AFE vehicle at 
its flight trajectory point at 77.8 km. The solutions 
were computed using an explicit, finite-rate, ther- 
mochemical nonequilbrium flow solver 12 . The same 
behavior as seen in the perfect gas computations is 
evident here. The cylindrical formulation does a far 
superior job in resolving flow quantities near the sin- 
gular line. 

Results Using the Roe-Upwind Solver 

Table 1. shows the various numerical experi- 
ments conducted to verify the accuracy of the Roe- 
upwind solver. As before, results from the invis- 
cid computations are compared with the results of 
Ref. 1. To simulate the axisymmetric flow around 
the spherical body, either a grid consisting of three 
planes separated by 1 deg. angle, shown in Fig. 9 
or a full three-dimensional grid around one-half of 
the body, similar to the grid shown in Fig. 1, is 
used. For inviscid cases, the grid points are placed 
uniformly between the body and the outer bound- 
ary and for viscous flows, the grid is stretched to- 
wards the wall. The singular line boundary exists in 
both the 3-plane grid and the full three-dimensional 
grid. Only explicit solutions were sought for this 
study and in all but one case, global time accuracy 


was maintained in obtaining solutions to the steady 
state. For the case where local time step was used, 
the step size in time was selected based on the maxi- 
mum local eigenvalue and the grid size such that the 
local Courant number was less than unity. For the 
viscous cases on a stretched grid, the flow near the 
body may restrict the time step size based on vis- 
cous stability criterion. However in the cases stud- 
ied here, the inviscid step size limitation from the 
shock region seems to determine overall time step 
size. Both the L^i^Q) und Max (A Q) norms were 
monitored during the iterative process and the con- 
vergence criteria of at least 4 orders of magnitude 
reduction in both these norms were required. The 
3-D code is vectorized and required 10 microseconds 
of CPU time on CRAY-YMP per grid point per it- 
eration. 

The computed solutions for the inviscid flow 
around the sphere geometry at a freestream Mach 
number of 20 are presented in terms of surface den- 
sity, pressure and shock standoff distances in Fig- 
ures 10 (a), (b) and (c) respectively. The coarse-grid 
first-order solution (case 1), coarse-grid second-order 
solution (case 2), fine-grid third-order accurate solu- 
tion (case 5) and the solution from Ref.l and are also 
compared in these Figures. The present computed 
solutions did not show any problem at or near the 
singular line region. The surface density, pressure 
and shock standoff distance comparisons between 
the present computations and that of Ref. 1 is found 
to be very good. The grid sensitivity analysis (case 
2 and 5) did not show any significant difference , 
due to sufficient second/third-order solution accu- 
racy, but the coarse-grid, first-order solution showed 
some differences, as seen in Figs. 10(a), (b) and (c). 

To compare the solution accuracy across the 
shock, normalized density (^), normalized pressure 
(-E-) and normalized u- velocity (~, where is 
the free-stream speed of sound) are shown in Figure 
11(a), (b), and (c) respectively. The first-order solu- 
tion, as expected, shows considerable shock spread 
compared to the higher-order solutions. The fine 
grid (case 5) solution captures the shock better and 
the shock location is also predicted well. The nor- 
malized density and pressure along 9 = 45° and 
9 = 75° lines were compared and are shown in 
Figure 12 (a), 12(b), 13(a) and 13(b) respectively. 
Though the overall comparison between the present 
computations and the solution from Ref. 1 is good, 
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the peak density and pressure behind the shock at 
0 = 75° is underpredicted by the present shock cap- 
turing computations relative to the shock fitting val- 
ues. This is to be expected. 

The use of local time step to achieve faster con- 
vergence to steady state solution is common tech- 
nique. The solution process, during convergence to 
steady state, is found to be non-conservative if lo- 
cal time step is employed. If the local maximum 
eigenvalues are significantly different (as the case is 
near the shock region, where max(A upj,treom ) /max 
^downstream ^ j s Q f or( j er G f f re estream Mach 

number) then the conservation is violated the most 
in such a region. Use of local time step may allow 
waves to travel at best possible speeds based on lo- 
cal conditions, but the solution is conservative and 
steady only if A Q at each and every grid point is 
guaranteed to be less than the convergence crite- 
rion. With the local time step method, moinitoring 
the local state variables at each grid point, such as 
the maxAQ norm, is necessary. In Figure 14, the 
stagnation density profiles from a apparently con- 
verged solution, based on £2 norm, computed with 
local time step is shown. The solution near the sin- 
gular line and near the shock has not converged com- 
pletely. Further convergence to steady state was 
found necessary and was achieved by reducing the 
maxAQ norm. Hence caution is advised in the use 
of local time step method. 

To further validate the present technique for 
viscous flows, an adiabatic and an isothermal wall 
boundary cases were computed with a viscous grid. 
The Reynolds number based on freestream condi- 
tions and the radius of the sphere was 3.7 x 10 6 and 
the wall temperature was set to freestream static 
temperature for the isothermal case. Surface pres- 
sure and shock standoff distances from the viscous 
computations are compared with the inviscid solu- 
tions in Figures 15 (a) and (b) respectively. The 
surface pressure has no significant effect due to the 
growth of the viscous boundary and shows smooth 
solution behavior near the singular line. The differ- 
ences observed in the shock standoff distance is due 
to the viscous layer at the wall. The isothermal wall 
conditions (with cold wall boundary) has the largest 
thermal and momentum boundary layer and shows 
the maximum deviation for the inviscid shock stand- 
off distance. In Figures 16 (a) and (b), density con- 
tours on the symmetry plane from the inviscid and 


isothermal cases are shown. The smooth variation 
of the contours near the singular line demonstrates 
the effectiveness of the present solution procedure. 

To further demonstrate that the singular line 
boundary does not present any problem, the 3-D, 
inviscid, Mach 20, flow around the sphere at zero 
and 15 deg. angles of attack were computed using 
the present solver. The grid used in both these cases 
was the same and was similar to the grid used in 
case (2) except 12 circemferential planes were con- 
structed to span one-half of the sphere. The singular 
line coincided with the freestream direction for the 
0 deg. angle of attack case and for the 15 deg. angle 
of attack case, the freestream and the singular line 
were 15 deg. apart. The intent here was not only to 
demonstrate the effectiveness of obtaining accurate 
solutions when there is flow across the singular line, 
but also to show that the solution symmetry regard- 
less of the location and orientation of the singular 
line. Since the geometry is a sphere, the 15 deg. an- 
gle of attack solution can be rotated by (-15) degrees 
and compared with the 0 deg. solution. 

In Figure 17(a) and (b), pressure contours on 
the bi-symmetry plane are shown for the 0 deg. and 
15 deg. angles of attack cases. In these figures, in 
addition to the pressure contour lines, the grid on 
the spherical surface and the singular line boundary 
are shown. The computed solution is very symm- 
teric about the stagnation point in both these cases. 
The surface pressure on the bi-symmtery plane are 
compared with the solution of Ref. 1 in Figures 18 
(a) and (b). The comparison is good and the lo- 
cation of the singular line has no effect on the so- 
lution. These two computations demonstrates the 
effectiveness of the present solution methodology to 
compute 3-D flows with grids that contain singular 
line boundaries. 

Concluding Remarks 

The methods presented in this study effectively 
reduce the difficulties associated with the singular 
line in 3-D blunt body flows. The first method refor- 
mulated the governing equations in terms of a mod- 
ified definition of the grid Jacobian. The resulting 
3-D solutions showed marked improvement over so- 
lutions when the usual, Cartesian grid Jacobian is 
used. The singularity is essentially moved to two 
source vectors which are not evaluated at the singu- 
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lar line. Both codes used identical grids, singular line 
boundary conditions, and first-order differencing of 
the inviscid fluxes. The Cartesian code underpre- 
dicted stagnation line shock standoff distance and 
placed the point of peak surface pressure away from 
the stagnation point. The cylindrical method was 
equally effective in improving viscous and real-gas 
solutions. 

The second method, which utilizes finite- 
volume techniques to properly evaluate of the met- 
ric and Jacobian of the transformation and accu- 
rate boundary condition procedures, demonstrated 
that the singular line problem can easily be avoided 
for axisymmteric and 3- D flow problems. Though 
the present solution procedure is found to work well 
with Roe’s upwind differencing technique and with 
MUSCL extrpolation with differential limiters, other 
forms of flux-splitting and flux-limiting schemes have 
yet to be explored. 

Both of the methods presented in this study il- 
lustrate that the key to effectively treating the singu- 
lar line boundary lies not in the boundary conditions 
themselves but rather in the proper determination of 
the metric and flux terms on the singular line. 
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Figure 8(a). Velocity magnitude contours 
- Cartesian formulation. 



Figure 8(b). Velocity magnitude contours 
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Case 

Grid Size 

Order 

Comments 

1 

25*3*25 

1st Order 

Inviscid 

2 

25*3*25 

2 & 3 Order 

Inviscid 

3 

49*3*49 

1st Order 

Inviscid 

4 

49*3*49 

2 & 3 Order 

Inviscid, variable/local time steps 

5 

49*3*49 

2 & 3 Order 

Inviscid. global time step 

6 

49*3*49 

3rd Order 

Viscous, adiabatic wall 

7 

49*3*49 

3rd Order 

Viscous, isothermal wall 

8 

25*12*25 

3rd Order 

Inviscid . a = 0.0 

9 

25*12*25 

3rd Order 

Inviscid, a= 15 

10 

25*3*25 

3rd Order 

Inviscid. adapted grid 


Table 1. Details of the cases run with the 3-D solver 
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Figure 9. Three-plane 3-D grid for 
axisymmteric flow 
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Figure 16(a). Density Contours on the 
symmetry plane , invscid conditions (case 5) 


Figure 16(b). Density Contours on the 
symmetry plane : viscous conditions (case 7) 


Figure 17(a). Pressure contours on the 

bi-symmtery plane - Case 8 


Figure 17(b). Pressure contours on the 
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Abstract 


Numerical simulations are used to investigate periodic combustion instabilities observed 
in ballistic-range experiments of blunt bodies flying at supersonic speeds through hydrogen- 
air mixtures. The computations are validated by comparing experimental shadowgraphs to 
shadowgraphs created from the computed flowflelds and by comparing the frequency of the 
experimental and computed instabilities. The numerical simulations utilize a logarithmic 
transformation of the species conservation equations which recently has been proposed as 
a way to reduce grid requirements for computing flows with shock-induced combustion. 
The transformation is applied to the Euler equations coupled to a detailed hydrogen- 
air chemical reaction mechanism with thirteen species and thirty-three reactions. The 
resulting differential equations axe solved using a finite- volume formulation and a two-step 
predictor-corrector scheme to advance the solution in time. Results are presented and 
compared for both a flux-vector splitting scheme and an upwind TVD scheme. The 
simulations add insight into the physical processes observed in the experiments and help 
establish the advantages of a logarithmic form of the species conservation equations for 
combustion flows. The usefulness of the ballistic-range experiments for the validation of 
numerical methods and chemical kinetic models is also demonstrated. 

Introduction 

Ballistic-range experiments performed in the 1960s and early 1970s 1-7 provide excellent 

data for studying the coupling between supersonic fluid dynamics and nonequilibrium 

\ 

chemical kinetics els well as for evaluating combustion flow codes. In these experiments, 
small projectiles were fired at supersonic speeds into a variety of premixed combustible 
mixtures. Shadowgraphs of the flowflelds exhibit two distinguishing features: one is the 
bow shock ahead of the projectile and the other is an energy-release front created by the 
ignition of the heated mixture behind the bow shock. Both features can be seen in the 
shadowgraph by Lehr 5 of a steady ballistic-range experiment at Mach 6.46. The region 
between the bow shock and the energy-release front is called the induction zone and it exists 
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because there is an ignition delay caused by chemical nonequilibrium. The induction zone 
is characterized by neax- constant values (the post-shock values) of temperature, density, 
pressure, and velocity. The size of the induction zone is determined by the fluid speed 
downstream of the shock and the ignition time corresponding to the post-shock conditions. 
When ignition occurs, the energy is released over an interval much shorter than the ignition 
delay time and appears as a discontinuity which is referred to as the energy- release front. 
Across the energy-release front the pressure is nearly constant, the temperature rises, and 
the density drops. Such behavior can be seen in the ballistic-range shadowgraph of Figure 1 
where the shift from light to dark across the bow shock corresponds to the density increase 
across the shock while the shift from dark to light across the energy-release front indicates 
a drop in density. 

Depending on the conditions of the experiment, one will observe either steady or 
unsteady flow. Projectile speeds greater than the detonation wave speed tend to induce 
steady flows while speeds less than the detonation wave speed can produce unsteady 
flows. The unsteady flows are characterized by two different regimes. One is called 
the regular regime and the other is called the large disturbance or irregular regime. All 
the unsteady simulations presented in this work are in the regular regime (see Alpert 
and Toong 7 for more on the large disturbance regime). Figures 2 a and 2b show an 
unsteady ballistic-range experiment with the same free stream conditions as in Figure 1 

(thg detonation wave speed is Mach. 5:11.) except that the projectile speed is Mach 4.79. 

t * ' % 

" \ 

The shadowgraphs reveal remarkable high-frequency oscillations. The frequency of the 

* s* 

observed instability is approximately 720 KHz as deduced from the the projectile speed 
and counting the number of oscillations which occur over a known distance. Figures 2a 
and 2b are from the same ballistic-range shot but show different view angles. The view 
axis of Figure 2a is perpendicular to the flight axis and Figure 2b is off-axis and reveals 
some of the three-dimensional structure of the flowfield. While being complex in many 
ways, the physics of these ballistic-range flows are predominantly driven by reaction 
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kinetics and convection phenomena; therefore, the complications and uncertainties of 
diffusion and mixing are removed from the problem. As a result, differences between 
the experimental data and numerical calculations can be attributed either to numerical 
errors or to improperly modeled chemical kinetics. 

An earlier work by Wilson and MacCormack 8 focused on the steady ballistic-range 
experiments and established the physical and numerical modeling requirements for accurate 
computations of these flows. It was shown that one of the most challenging aspects 
is providing adequate resolution in the induction zone between the bow shock and the 
energy-release front. Fine grid resolution is required in this region not because there are 
large changes in the flow quantities such as pressure, density, or velocity (these quantities 
are nearly constant in the induction zone) but because the mass-fractions of the important 
radical species change exponentially in this region. Adaptive grid techniques offer one 
way to efficiently distribute points so that grid resolution is provided where it is needed 
most. These techniques, however, are not always adequate. The induction zone may 
cover such a large region of the flowfield that grid requirements still remain large. Also, 
complex flow structures such as those seen in the unsteady ballistic-range cases require 
very sophisticated (and probably very expensive) grid adaption techniques. Furthermore, 
the small grid spacings in the adapted regions result in undesired time step limitations. 

To avoid small grid spacing in the induction zone, Sussman and Wilson 9 proposed a 
logarithmic transformation of the species conservation equations. The reasoning for this 

% • ’ S 

transformation is discussed below. Figure '3'shows the mass fraction of the hydroxyl radical 

• N 

(OH) in the induction zone for the case of shock-induced combustion which was presented 
by Sussman and Wilson. 9 The exponential growth of the hydroxyl mass fraction is seen as 
linear growth on the logarithmic scale. Many grid points are required to accurately capture 
this growth because a finite-difference approximation for a derivative of an exponentially 
varying function can result in large errors. In contrast, fewer points are needed to resolve 
the linear behavior seen on the logarithmic scale. 
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The use of the logarithmic transformation allows accurate simulations of both the steady 
and unsteady ballistic-range cases on grids more coarse than have been used in previous 
works. The reduced grid densities make it possible for the first time to simulate the 
unsteady ballistic-range cases with a detailed chemical kinetic model. The chemical kinetic 
model used in this work can be found in Wilson and MacCormack. 8 It is the mechanism 
proposed by Jachimowski 10 containing 13 chemical species and 33 reactions except that 
one of the reaction rate expressions (the expression most influential for ignition times) has 
been replaced by an expression recommended by another source. This modification was 
found by Wilson and MacCormack 8 to give better agreement with experiment. Additional 
details on the chemical kinetics are found in Wilson. 11 

The simulations of the unsteady flows in this work add insight into the combustion insta- 
bilities observed in ballistic-range experiments and allow previously proposed mechanisms 
for the instabilities to be assessed. In particular, this work considers the wave interaction 
mechanism of McVey and Toong 6 and the modified mechanism proposed by Matsuo and 
Fujiwara 12 . McVey and Toong developed their mechanism by using experimental data and 
analytical calculations. Matsuo and Fujiwaxa used modern computational fluid dynamics 
techniques, fine grids, and a two-equation, two-step combustion model to do qualitative 
simulations which led them to propose their modified mechanism. The use of detailed 
combustion kinetics in this work allows quantitative studies and direct comparison to 

experiment to be made. The experimentally measured oscillation frequencies are compared 

\ 

to’ the values predicted by numerical simulations. In addition to adding greater physical 
understanding, the current method also provides a tool to validate and possibly “tune” 
proposed chemical kinetic models. 

Logarithmic Form of the Species Conservation Equation 

This section gives an overview of the logarithmic transformation of the species conser- 
vation equations as presented by Sussman and Wilson. 9 The transformation is based on 
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the definition 


7r s = pln(y) = /oln(c a ), (1) 

where p is the total density, p s is the density of species 5 , and c s is the mass fraction of 
species s. When c 3 is a exponentially varying quantity, tt 3 is a linearly varying quantity 
(assuming constant p). A special property of the relationship in Equation (1) is that it 
can be used along with the total mass conservation equation, 

Jt p + £] pUi = 0 ’ (2) 

to transform the species mass conservation equation (neglecting diffusion), 

d d 

Wt fi, + - (3) 

into a transformed species equation, 


d d 

-X-TI-s + -r TT 3 Uj — w , 

at oxj 


P_ 

Ps’’ 


( 4 ) 


where Uj is the j th component of velocity and w 3 is the chemical source term for species 
s. The transformation from Equation (2) to (4) retains conservation form. 


In theory, all of the species mass conservation equations can be replaced with Equa- 
tion (4). In practice, however, this is not done. Although Equation (4) is in conservation 
form, it does not guarantee conservation of mass or conservation of the individual elements. 
Enforcement of these constraints is achieved by replacing some of the logarithmic equations 
with' element mass conservation equation^ which are written in the formi 

d * v d * 

Jt p ‘ + d7j P ‘ Ui = °’ (5) 

where p * is the total'density of element e contained in all of the species (the * is used to 
avoid confusion between element and species densities; p* N ^ Pn)- There are as many of 
these equations as there are elements in the chemical model. The total density is given by 

lem&nta 

P= E rt- (6) 

e=l 
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The chemical reaction mechanism used in this work involves 13 species: JV 2 , O 2 , H 2 , 
NO, OH, NO 2 , HNO, HO 2 , H 2 O, H 2 O 2 , N, O, H. Since three elements are present 
in the mechanism (iV, 0, H ), three equations of the form of Equation (5) are required. 
The densities p* N , p~Q, and p* H are defined using expressions which are the sum of the 
contributions from each of the species containing the element of interest. For example, the 
expression for p* N depends on JV 2 , NO, NO 2 , HNO , and N and is written 

Pj± _ o PN? PNO PNO2 PHNO , _PN_ /-x 

Mat Mn 2 Mno Mho 2 Mhno Mn 


where M s is the molecular mass of species s. The choice of which three species conservation 
equations to replace with element mass conservation equations is not unique.- The selection 
of N 2 , O 2 , and H 2 seems appropriate because these species d"o not exhibit exponential 
growth of mass fraction. The removal of these 3 species conservation equations means that 
the corresponding 3 species densities do not appear in the set of conserved variables. They 
can, however, be written in terms of the conservation variables with the equations of the 
type in Equation (7). The remaining ten species conservation equations are written in the 
logarithmic form of the mass conservation equation, Equation (4) (for species NO, OH, 
N0 2 , HNO, H0 2 , H 2 0, H 2 O 2 , N, 0, and H). 


Numerical Formulation 


This section presents the Euler equations . with the logarithmic form of the species 

, x ' 

conservation equations. The 13 species ,mass conservation equations of the conventional 
formulation have been replaced by-equations of similar form except that 3 are elemental 
mass conservation equations (See Equation (5)) and 10 are in the logarithmic variable tt s 
(See Equation (4)). The momentum and energy equations in the logarithmic formulation 
are the same as they were in the conventional formulation. In vector notation the governing 
equations are written 


dU_ dF_ dG __ 

dt dx dy 


( 8 ) 
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where U is the state vector of conserved quantities, F and G axe the convective (inviscid) 
flux vectors in the x and y coordinate directions, respectively, and W is the vector of source 
terms. These vectors axe given below: 

P*N \ / Pn u \ / P*n v \ 0 

Po Po u Po v I 0 

Ph Ph u P*h v 0 

7T/VO KNOU KNOV wno/cno 

' , F= 1 , G= 1 , W= [ 

TTH KHU KhV WH/CH 

pu pu 2 + p puv 0 

pv puv pv 2 + p • o - 

E v uE v vE v __ w v 

E / \u(E + p )) \v(E + p)J \ 0 

where E v is vibrational energy per unit volume, w v is the- source term representing 
vibrational relaxation, E is the total energy per unit volume, and p is pressure. A separate 
vibrational energy equation is not required for the present computations but has been 
retained from previous work. For the cases presented here, thermal equilibrium is simulated 
by using a small time constant for vibrational relaxation. 

In Equation (8), pressure is a homogeneous function of degree one in the conserved 
variables U and the flux Jacobians (| jj and have the same eigenvalues as that of the 
formulation with the conventional form of the species mass conservation equations. There- 
fore, existing numerical techniques can be applied in a straightforward manner. In this 
work, two commonly used shock-capturing numerical techniques axe employed. One is a 
modified Steger-Waxming flux-vector splitting technique based on work by MacCormack 13 
and Candler 14 (this method shall be referred to as the flux-vector splitting scheme) and the 
other is the Harten-Yee upwind TVD (non-MUSCL) scheme 15 (referred to as the upwind 
TVD scheme). Both techniques axe spatially second order and solved in a point-implicit 
manner. That is, the convective terms axe treated explicitly whereas the source terms are 
treated implicitly in time (because of the small time scales of the combustion chemistry 
and vibrational relaxation). Therefore, a block inversion is required at each point at each 
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time step. Details on the flux-vector splitting scheme using the logarithmic formulation 
can be found in Wilson. 11 

Numerical Simulations 

This section presents numerical simulations of ballistic-range shadowgraphs by Lehr. 4 ’ 5 
The particular experiments of interest used spherical-nosed projectiles with cylindrical 
afterbodies of 15 millimeters diameter. The cases include a range of Mach numbers so 
that both steady and unsteady flows are represented. All the cases to be considered 
have freestream with a temperature of 292 K, a pressure of 320 mm Hg, and a premixed 
stoichiometric hydrogen-air mixture. At these conditions, the detonation wave speed is 
Mach 5.11. 

Before unsteady simulations are presented, a simulation of the steady Mach 6.46 case 
in Figure 1 is shown to demonstrate the advantage of the logarithmic formulation and to 
give an understanding of why it is practical to use this formulation for the unsteady cases. 
Figure 4 contains density contours for computations of the Mach 6.46 case using flux- vector 
splitting and the upwind TVD methods, respectively. Both computations use the 52 x 52 
grid shown in Figure 5. The experimental bow shock and energy- release front positions are 
overplotted on the density contours and show that both numerical methods provide good 
agreement with the experiment. A numerical simulation using the conventional form of the 
species conservation equations can be found in Wilson and MacCormack. 8 Good agreement 

with experiment was found in that work as well but a 321 x 64 (321 points around the 

* * \ * 

body and 64 points normal to the body) 'mesh with grid adaption was required. The grid 
for the calculation using the logarithmic formulation is more than 8 times smaller because 
the induction zone is- adequately resolved with fewer grid points. 

Simulations of Lehr’s unsteady Mach 4.79 case shown in Figure 2 are now presented. The 
computations were started from a flowfield that was initially set to freestream conditions 
everywhere and advanced in time without chemistry until a bow shock was established 
at which time the combustion chemistry was turned on. No other special procedures 
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were used. The computations used a nearly equally-spaced 375 x 161 grid. The fine grid 
was not needed to capture the unsteady behavior but was used to resolve some of the 
finer structures of the flow. Computation on a 375 x 81 grid yielded similar overall flow 
features and an oscillation frequency near that predicted by the 375 x 161 grid. 11 Density 
contour plots for both the flux-vector splitting and upwind TVD methods are presented in 
Figures 6a and 66, respectively. The figures represent one point in time and one point in the 
period of the instability (the two solutions are not at exactly corresponding points within a 
period) and show that both computational methods predict unsteady behavior. As in the 
steady computation of Figure 4, the density contours exhibit an outer bow shock followed 
by an induction zone and an energy- release front. The most obvious difference between 
the unsteady calculation and the steady one is the pulsing structure of the energy- release 
front. These pulses are seen in the experimental shadowgraph with the off-axis view in 
Figure 26. 

Figures 7 a and 76 contain shadowgraphs computed from the flux-vector splitting and 
upwind TVD flowfields, respectively. 16 A comparison of these computed shadowgraphs 
with the enlargement of Figure 2a presented in Figure 8 clearly shows many similarities. 
Note that the pulses in the energy-release front create the vertical line pattern when the 
axisymmetric flowfield is projected onto a plane. There is a larger streamwise separation 
between the vertical lines in the computed shadowgraph compared to the experiment and 
thus the computation predicts a lower frequency than the experiment.. Both numerical 
methods predict an instability, frequency' of approximately 530 KHz compared to the 

* s' 

measured value of 720 KHz. The primary difference between the flux-vector splitting 
and upwind TVD solutions is the resolution of the flow structures near the shoulder region 
of the projectile. The flux-vector splitting smears the features in this region to a larger 
degree. 

Further validation of the computations is provided by comparing the smaller features 
in the flow. A schematic of the flowfield which labels the different structures is found 
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in Figure 9. A thorough analysis of these structures by McVey and Toong 6 identified 
the wave extending from the maximum amplitude of each energy-release front pulsation 
to the bow shock as a contact discontinuity (McVey and Toong refer to it as an entropy 
wave). Its identity was established by observing that the wave remained at a fixed position 
on a pulsation as the pulsation was convected downstream. A compression or expansion 
wave would have moved faster than the local flow with a relative speed which depends on 
the speed of sound in the mixture. The wave that looks like a reflection of the contact 
discontinuity off the bow shock does move at a speed greater than the . local flow and 
is identified as a compression wave. The dark-light shadings across the waves further 
support the asserted identities of the waves. All of these features can be observed in both 
the experimental and simulated shadowgraphs. 

Perhaps the most convincing observation supporting the accuracy of the numerical 
simulation is the presence of some of the weakest observable flow structures in both 
the experiment and computation. Although the literature only mentions one contact 
discontinuity originating from the tip of each pulsation, two contact discontinuities can 
be observed if the flowfield is examined closely. The additional contact discontinuity is 
weaker than the first and is at a more shallow angle. Both of these features are observed 
in the computed shadowgraph as well as in the experimental shadowgraph. 

Investigation of other flowfield quantities adds information that cannot be deduced 
from, the shadowgraphs. Pressure contours reveal that each pulse of the energy- release 

N * 

front creates a compression wave which 'travels both upstream toward the bow shock and 
downstream toward the projectile body. Figure 10 shows the pressure contours at one 
point in time. The wave patterns between the bow shock and the projectile appear quite 
complicated but are simplified by noting that there are primarily two families. One family 
is made up of compression waves originating directly from the energy-release front and the 
other family is made up of compression waves from the first family which have reflected 
off the projectile body. 


11 



Instability Mechanism 


The phenomenon which causes and sustains the combustion instabilities observed in 
the ballistic-range has not been directly observed. Therefore, plausible explanations for 
the unsteadiness have been developed by extrapolating data from outside the nose region 
(mostly from the shadowgraphs) and by transferring knowledge from other flow situations. 
It was not until the work of Toong and his associates that plausible detailed mechanisms 
for the regular 6 and large-disturbance 7 regimes were put forward. 

The work of McVey and Toong 6 is discussed further here because it is supported by the 
present numerical simulations. McVey and Toong used what they called a wave-interaction 
model to explain the oscillations observed in the regular regime. All of .the important 
processes in this model occur between the energy-release front and the bow shock front in 
the stagnation region of the nose. As an aid to understanding the model, several of the 
primary components of the model are isolated and described first. 

A major component in the McVey and Toong model is the interaction between the bow 
shock and a compression wave. Figure 11 depicts a one-dimensional flow at several different 
times (time proceeds from bottom to top). Initially there is a stationary normal shock and 
a compression wave downstream of the shock moving toward the shock. At a later time 
the compression wave overtakes the shock, strengthening it and causing the shock to move 
forward. The strengthened shock causes a change in the fluid properties behind it. Most 
important for the instability is a higher fluid temperature. This higher temperature fluid is 
separated from the fluid which has crossed the shock at the original strength by a contact 

- s' 

discontinuity which convects downstream at the fluid velocity. Additionally, a rarefaction 
wave is created which-travels downstream. 

Interesting, flow features are created when the interaction of Figure 11 is combined with 
the flow of a combustible mixture. This is the situation presented in Figure 12a. As in 
Figure 11, there is initially a stationary normal shock wave with a compression moving 
toward it from the downstream side. Now, there is also burned fluid downstream of the 


12 



shock at some induction length. As before, the compression wave overtakes the bow shock, 
strengthening it and creating a contact discontinuity and a rarefaction wave. Because the 
fluid on the upstream side of the contact discontinuity is hotter it has a shorter induction 
time and burns more quickly. The result is that for a time there are two regions of 
burning, one at an induction length corresponding to the conditions behind the original 
shock strength and one determined by the shorter induction length of the hotter fluid 
behind the strengthened bow shock (McVey and Toong show that the effect of the contact 
discontinuity is more important than the effect of the compression wave or rarefaction wave 
on the burning locations). Figure 126 is an x-t diagram corresponding to the interaction 
in Figure 12a. The general features of this diagram are recognizable in the McVey and 
Toong mechanism presented next. 

The McVey and Toong mechanism is schematically depicted in the x-t diagram in 
Figure 13. The diagram contains the features along the stagnation streamline in time 
and is explained in the following four steps: 

1) A compression wave on the downstream side of the bow shock approaches the bow 
shock. The compression wave overtakes the bow shock causing the bow shock to 
move forward, thus creating a reflected rarefaction (which is weak and is ignored in 
the model) and a contact discontinuity. The gas on the upstream side of the contact 
discontinuity is hotter than on the downstream side because it has passed through the 
•strengthened bow shock. ■ \ • 

/ * 1 * N ' 
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. 2) The rarefaction wave propagating upstream from the energy-release front overtakes 
the bow shock, wakening it and restoring it to its original strength. The origin of the 
rarefaction wave is discussed in Step 4. 

3) The contact discontinuity created in Step 1 convects downstream. At a point between 
the bow shock and the energy-release front, a new energy- release front is produced 
because the hotter gas on the upstream side of the contact discontinuity is ignited more 
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quickly than the gas on the downstream side of the discontinuity. Because the ignition 
delay time has an exponential (Arrhenius) dependence on temperature, the change 
in ignition time is significant even with small temperature changes across the contact 
discontinuity. The creation of a new zone of combustion, in turn, creates compression 
waves which propagate both upstream towaxd the bow shock (this is the wave which 
overtakes the bow shock in step 1) and downstream toward the projectile. It can be 
shown that the upstream and downstream compression waves are necessary to match 
the fluid velocity jumps across the new energy- release front and satisfy conservation 
of momentum (see Alpert and Toong 7 ). 

4) The contact discontinuity (with burning on the upstream side) eventually reaches the 
position of the original energy- release front. This extinguishes the original energy- 
release front and creates rarefaction waves which propagate in the upstream and 
downstream direction. The energy-release front then begins to recede toward the 
location of the original front as colder gas (through a the bow shock of original strength) 
reaches the front. 

McVey and Toong describe how the compression waves and contact discontinuities present 
in the interaction model can explain the various features observed in experimental shad- 
owgraphs. From Figure 13 it is seen that the period of an oscillation is the time required 
for the contact discontinuity to travel from the shock to the reaction front and then for 
an acoustic wave to travel back to the bow shock. If the post-shock conditions are known, 
an analytical expression for the period of. observed experimental oscillations can be found. 
Given a period of oscillation provided from an experiment, McVey and Toong were able 
to use the expression, to predict the induction time for the gas mixture. Their predicted 
induction times agreed with previously published data. 

Figure 14 contains computed x-t diagrams of density and pressure along the stagnation 
streamline for the Mach 4.79 case from the computations presented earlier. A comparison 
of these computed diagrams with the McVey and Toong mechanism in Figure 13 reveal 
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many similarities. The density contours show a changing position of energy-release front in 
time that is identical to the pattern predicted by the wave interaction model (i.e. note that 
the jagged edge of the energy-release front appears in both Figures 13 and 14a). In the 
region between the bow shock and the energy-release front, density contours also show that 
a contact discontinuity is created at the bow shock when a compression wave originating 
from the energy-release front overtakes it. The hotter gas on the upstream side of the 
contact discontinuity reduces the ignition delay time and causes a new energy-release front 
to be created. The new energy-release front is the source of the compression waves traveling 
in both the upstream and downstream direction. These compression waves are clearly seen 
in the pressure contours of Figure 146. The identity of the contact discontinuities labeled 
in Figure 14a is verified by their absence in the pressure contours of Figure 146 because a 
contact discontinuity has no pressure jump across it. 

A phenomenon seen in the computations which was apparently not anticipated by 
McVey and Toong is the existence and/or importance of the compression waves reflecting 
off the projectile nose. Figure 146 shows that a compression wave created at the energy- 
release front travels toward the projectile nose, reflects, and eventually overtakes the 
bow shock. In this particular calculation, the compression wave reflects off the nose 
and moves toward the bow shock overtaking it at nearly the same time that the most 
recently created compression wave arrives at the bow shock. The coordination between the 

compression wave and the reflected compression wave is not always exact and therefore the 
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two waves sometimes hit the bow shock, at slightly different times. This tends to slightly 
smear or deform some of the structures. The importance of considering the reflection 
of the compression waves off the projectile has been also been recognized by the recent 
calculations of Matsuo and Fujiwara. 12 

Figure 15 contains computed x-t diagrams of density and pressure from a case at 
Mach 5.04. This case is interesting because the interaction of the wave reflected off the 
projectile is significantly different than is seen in the Mach 4.79 case. As in the Mach 4.79 
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case, the interactions occuring between the bow shock and energy-release front corresponds 
to the model proposed by McVey and Toong. The contact discontinuities are clearly defined 
by the density contours and are absent from the pressure contours. This case differs from 
the Mach 4.79 case because during the time that a compression wave travels from the 
energy release front to the projectile body and back, two new compression waves axe 
created (i.e. two periods have passed). The computed frequency is 820 Khz while the 
measured frequency is 1.04 MHz. A thorough discussion of this case is given in Wilson 11 
along with a comparison to an experiment at Mach 4.18 with a much lower oscillation 
frequency. 

Frequency Sensitivity to the Hydrogen-Air Reaction Mechanism 

The nearly constant values of the instability frequency predicted by the numerical 
simulations using different grid sizes (375 x 321, 375 x 161, and 375 x 81) and two 
different numerical schemes suggests that the underprediction of the frequency in the 
computations is not caused by numerical error but by the chemical kinetic model. Wilson 
and MacCormack 8 established that the blunt body exothermic flow fields are quite sensitive 
to the chain- branching reaction 


H + O2-+OH + O. (Rl) 

In that work, flowfields were computed with two different rate expressions for this reaction. 
One-was the original expression from the Jachimowski 10 reaction model and the other was 
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from Warnatz 17 . The rate expression .recommended by Warnatz was adopted because it 
gave better agreement with experiment for the steady cases (this is the rate expression 
that has been used for all cases presented in this paper). To investigate the influence 
of the chemical kinetic mechanism on the frequency of the oscillations, the Mach 4.79 
simulation is repeated using the Jachimowski rate expression for Reaction (Rl). This 
rate expression gives shorter ignition times at high temperatures (T > 1400/v) than 
the rate constant recommended by Warnatz. This change should lead to an increase 
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in the oscillation frequency because the induction zone is smaller and therefore the travel 
time between the energy-release front and the bow shock is reduced. Consistent with 
this expectation, the computed frequency of the oscillations increased from approximately 
530 KHz to 820 KHz. Since the experimentally determined frequency is 720 KHz, it is 
concluded that the uncertainties in the rate constants for the reaction mechanism could 
explain the differences between the experiments and the computations. The sensitivity of 
the oscillation frequency to the chemical reactions make the simulations a useful validation 
tool for kinetic models. The unsteady cases appear to provide a better measure of a kinetic 
model than the steady cases because with the available data the oscillation frequency can 
be determined with greater precision than the positions of the bow shock and energy- release 
front. 

As a final note, it is mentioned that the arrival of new compression waves and reflected 
pressure waves at the bow shock sometimes differ significantly (this occurred in the Mach 
4.79 case with the Jachimowski rate expression 11 ). The separate interactions with the bow 
shock by each wave lead to separate contact discontinuities and cause the x-t diagrams 
to appear less organized. It seems, however, that even with these waves out of phase, 
an overall periodic nature can exist because the reflected waves have a secondary effect. 
Further investigations are required to assess the precise contribution to the combustion 
instability by the compression waves reflected off the projectile nose. It is possible that the 
reflected waves are more important to the instabilities seen in the large-disturbance regime 
investigated by Alpert and Toong 8 ’ than 'they are to the regular regime cases considered 
here. 

Conclusion 

The successful simulations of both steady and unsteady exothermic ballistic-range 
experiments have added to the understanding of high speed flows with combustion and have 
led to improvements in the numerical techniques for these types of flows. Namely, proper 
grid resolution is a primary challenge of simulations with combustion phenomena and these 
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grid requirements axe reduced significantly by the use of the logarithmic transformation of 
the species conservation equations. The simulations have also provided information about 
combustion instabilities which may, in turn, help interpret unstable combustion phenomena 
in other settings. The accuracy of the simulations is supported by grid refinement studies 
and the use of two different numerical methods, a flux vector splitting scheme and an 
upwind TVD scheme. Further validation is provided by comparison to past analytical and 
numerical predictions of the instabilities. Finally, the simulations have reaffirmed the value 
of ballistic-range experiments as a source of data for high speed flows with combustion and 
have demonstrated their usefulness for validating numerical methods and chemical kinetic 
mechanisms. It was shown that differences in the experimental and computed instability 
frequencies are probably due to uncertainties in the chemical reaction rates and therefore 
the experiments provide a means of validating proposed chemical kinetic models. 
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Captions 


Figure 1: Shadowgraph of a Spherical Nose Projectile Moving at Mach 6.46 in a Stoichio- 
metric Hydrogen- Air Mixture (Courtesy of H.F. Lehr). 

Figure 2: Shadowgraph of a Spherical Nose Projectile Moving at Mach 4.79 in a Stoichio- 
metric Hydrogen- Air Mixture: a) side view b) off-axis view (Courtesy of H.F. Lehr). 

Figure 3: Variation of the mass fraction of the hydroxyl (OH) radical in the induction 
zone for a quasi-one-dimensional calculation of shock-induced combustion. • 

Figure 4: Density contours from a numerical simulation of the Mach 6.46 experiment in 
Figure 1: a) flux-vector splitting b) upwind TVD. 

Figure 5: 52 x 52 grid used for the calculations presented in Figure 4. 

Figure 6: Density contours from a numerical simulation of the Mach 4.79 experiment in 
Figure 2: a) flux- vector splitting b) upwind TVD. 

Figure 7: Computed shadowgraphs from a numerical simulation of the Mach 4.79 exper- 
iment in Figure 2: a) flux-vector splitting b) upwind TVD (shadowgraphs computed by 
L.A. Yates* 6 ). 

Figure 8: Enlargement of the M a c h 4 . 7 9_ b a 1 1 i s t i c - r a n g e shadowgraph in Figure 2a (cour- 
teSy of H.F. Lehr). ' ' 

Figure 9: Schematic; of the periodic flowfield structures seen in shadowgraphs of ballistic- 
range experiments. 

Figure 10: Pressure contours for the Mach 4.79 case at one point in an oscillation period. 

Figure 11: Schematic of a one-dimensional bow shock - compression wave interation. 
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Figure 12: Schematic of a bow shock - compression wave interaction with shock-induced 
combustion: a) series of one-dimensional diagrams b) x-t diagram. 

Figure 13: McVey and Toong 6 wave interaction model. 

Figure 14: Computed x-t diagrams of density and pressure along the stagnation streamline 
for the Mach 4.79 case using flux-vector splitting: a) density contours b) pressure contours. 

Figure 15: Computed x-t diagrams of density and pressure along the stagnation streamline 
for the Mach 5.04 case using flux-vector splitting: a) density contours b) pressure contours. 
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